Carrega de les llibreries que es necessiten
library(arules)
library(C50)
library(caret)
library(class)
library(cluster)
library(dplyr)
library(factoextra)
library(FactoMineR)
library(fpc)
library(ggpubr)
library(ggplot2)
library(grid)
library(gridExtra)
library(mclust)
library(plyr)
library(purrr)
library(randomForest)
library(rattle)
library(RColorBrewer)
library(rpart.plot)
library(rpart)
library(vcd)Aquesta pràctica cobreix de forma transversal l’assignatura.
Les Pràctiques 1 i 2 de l’assignatura es plantegen d’una forma conjunta de manera que la Pràctica 2 serà continuació de la 1.
L’objectiu global de les dues pràctiques consisteix en seleccionar un o diversos jocs de dades, realitzar les tasques de preparació i anàlisi exploratòria amb l’objectiu de disposar de dades llestos per aplicar algoritmes de clustering, associació i classificació.
A partir de la pràctica 1 vam generar des de 2 datasets originals (retail i bank) 5 datasets resultants que farems servir per als diferents exercicis de la pràctica 2.
bank-clean-discret: conjunt de dades obtingut a partir del dataset de marqueting bancari on s’han estat discretitzats alguns camps, eliminat alguns camps nulls i eliminar outliers.
bank-discret: conjunt de dades obtingut a partir del dataset de marqueting bancari on s’han estat discretitzats alguns camps. No s’han eliminat nulls ni outliers.
retail-customer: conjunt de dades obtingut del dataset de venta online on s’han generat nous atribtuts per a tenir dades de comportament dels clients.
retail-product: conjunt de dades obtingut del dataset de venta online on s’han generat nous atribtuts per a tenir dades dels productes.
retail-association: conjunt de dades obtingut del dataset de venta online on s’han generat transaccions de les compras realitzades.
A partir del dataset que hem preparat en la pràctica 1 anem a generar regles d’associació.
retail.association <- read.csv('retail-association.csv', sep = ',')
summary(retail.association)## InvoiceNo StockCode Description
## Min. :536365 85123A : 1686 WHITE HANGING HEART T-LIGHT HOLDER: 1680
## 1st Qu.:549316 85099B : 1329 JUMBO BAG RED RETROSPOT : 1329
## Median :562213 47566 : 1275 PARTY BUNTING : 1275
## Mean :560833 20725 : 1207 LUNCH BAG RED RETROSPOT : 1207
## 3rd Qu.:572302 84879 : 1159 ASSORTED COLOUR BIRD ORNAMENT : 1159
## Max. :581587 22720 : 1120 SET OF 3 CAKE TINS PANTRY DESIGN : 1120
## (Other):330406 (Other) :330412
## Quantity InvoiceDate UnitPrice CustomerID
## Min. : 1.000 11/14/2011 15:27: 460 Min. :0.000 Min. :12347
## 1st Qu.: 2.000 11/28/2011 15:54: 450 1st Qu.:1.250 1st Qu.:13994
## Median : 6.000 12/5/2011 17:17 : 442 Median :1.650 Median :15245
## Mean : 7.477 10/31/2011 14:09: 384 Mean :2.192 Mean :15326
## 3rd Qu.:12.000 11/23/2011 13:39: 370 3rd Qu.:2.950 3rd Qu.:16818
## Max. :27.000 9/21/2011 14:40 : 370 Max. :7.500 Max. :18287
## (Other) :335706
## Country Cancellation TotalPrice
## United Kingdom:305151 Mode :logical Min. : 0.00
## Germany : 7465 FALSE:338182 1st Qu.: 3.80
## France : 6905 Median : 10.08
## EIRE : 5451 Mean : 12.77
## Spain : 2046 3rd Qu.: 17.40
## Belgium : 1660 Max. :178.80
## (Other) : 9504Fem una transformació per obtenir les transaccions de les compres segons el InvoiceNo. El camps que ens interessa veure és Description.
transations_split <- split(x=retail.association[,"Description"],f=retail.association$InvoiceNo)
# trasformem el data com a transactions. Cada transacció es una compra (InvioceNo) i els productes que ha comprat
transations <- as(transations_split,"transactions")## Warning in asMethod(object): removing duplicated items in transactionsitemFrequencyPlot(transations, support = 0.04, cex.names = .6, col = rainbow(15), topN=30)
Mostrem els productes que més es repeteixen. En concret estem mostrant aquells que tenen un suport o freqüència >= 0.04 Això vold ir que, si tenim 16835 transaccions, estem mostrant x >= 16835 * 0.04, és a dir aquells que apareixen com a mínim 673,4 vegades.
Si llancem l’algoritme “apriori”, generarem directament un set de regles amb diferent support, confidence i lift.
El support indica quantes vegades s’han trovat les regles {lsh => rhs} en el dataset, com més alt millor. Es la “popularitat” d’un conjunt d’elements del dataset. On {lsh => rhs} indica que si es compra lsh, compra rhs.
La confidence ens parla de la probabilitat de comprar {rhs} si es compra {lhs} (lhs => rhs / lhs).
El lift és un paràmetre que ens indica quan d’aleatorietat hi ha a les regles. Un lift d’1 o menys ens indica que la regla és completament fruit de l’atzar. El lift ens diu eb una regla, com es produeix la regla en funció dels elements de la regla (support(lhs => rhs) / support(lhs) * support(rhs)).
Generem les regles per un support mínim de 0.02, confidence mínim de 0.4
rules <- apriori(transations, parameter = list(support = 0.02, confidence = 0.4))## Apriori
##
## Parameter specification:
## confidence minval smax arem aval originalSupport maxtime support minlen
## 0.4 0.1 1 none FALSE TRUE 5 0.02 1
## maxlen target ext
## 10 rules FALSE
##
## Algorithmic control:
## filter tree heap memopt load sort verbose
## 0.1 TRUE TRUE FALSE TRUE 2 TRUE
##
## Absolute minimum support count: 336
##
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[3578 item(s), 16835 transaction(s)] done [0.21s].
## sorting and recoding items ... [198 item(s)] done [0.00s].
## creating transaction tree ... done [0.01s].
## checking subsets of size 1 2 3 done [0.00s].
## writing ... [51 rule(s)] done [0.00s].
## creating S4 object ... done [0.00s].inspect(head(sort(rules, by = "confidence"), 10))## lhs rhs support confidence lift count
## [1] {PINK REGENCY TEACUP AND SAUCER,
## ROSES REGENCY TEACUP AND SAUCER } => {GREEN REGENCY TEACUP AND SAUCER} 0.02197802 0.8915663 22.40227 370
## [2] {GREEN REGENCY TEACUP AND SAUCER,
## PINK REGENCY TEACUP AND SAUCER} => {ROSES REGENCY TEACUP AND SAUCER } 0.02197802 0.8371041 18.84044 370
## [3] {PINK REGENCY TEACUP AND SAUCER} => {GREEN REGENCY TEACUP AND SAUCER} 0.02625483 0.8215613 20.64326 442
## [4] {GREEN REGENCY TEACUP AND SAUCER} => {ROSES REGENCY TEACUP AND SAUCER } 0.03070983 0.7716418 17.36710 517
## [5] {PINK REGENCY TEACUP AND SAUCER} => {ROSES REGENCY TEACUP AND SAUCER } 0.02465102 0.7713755 17.36110 415
## [6] {GARDENERS KNEELING PAD CUP OF TEA } => {GARDENERS KNEELING PAD KEEP CALM } 0.02583903 0.7213930 17.12927 435
## [7] {GREEN REGENCY TEACUP AND SAUCER,
## ROSES REGENCY TEACUP AND SAUCER } => {PINK REGENCY TEACUP AND SAUCER} 0.02197802 0.7156673 22.39453 370
## [8] {ROSES REGENCY TEACUP AND SAUCER } => {GREEN REGENCY TEACUP AND SAUCER} 0.03070983 0.6911765 17.36710 517
## [9] {DOLLY GIRL LUNCH BOX} => {SPACEBOY LUNCH BOX } 0.02257202 0.6713781 17.82752 380
## [10] {ALARM CLOCK BAKELIKE GREEN} => {ALARM CLOCK BAKELIKE RED } 0.03059103 0.6679637 13.21406 515Amb aquesta ordenació per confidence, tenim una probabilitat alta que si compren {PINK REGENCY TEACUP AND SAUCER, ROSES REGENCY TEACUP AND SAUCER} compraran {GREEN REGENCY TEACUP AND SAUCER}. És un exemple, el mateix passa amb la resta de les 9 regles mostrades. En tots els casos el lift és molt alt, per tant no hi ha cap tipus d’aleatorietat en els resultats.
inspect(head(sort(rules, by = "support"), 10))## lhs rhs support confidence lift count
## [1] {GREEN REGENCY TEACUP AND SAUCER} => {ROSES REGENCY TEACUP AND SAUCER } 0.03070983 0.7716418 17.367098 517
## [2] {ROSES REGENCY TEACUP AND SAUCER } => {GREEN REGENCY TEACUP AND SAUCER} 0.03070983 0.6911765 17.367098 517
## [3] {ALARM CLOCK BAKELIKE GREEN} => {ALARM CLOCK BAKELIKE RED } 0.03059103 0.6679637 13.214064 515
## [4] {ALARM CLOCK BAKELIKE RED } => {ALARM CLOCK BAKELIKE GREEN} 0.03059103 0.6051704 13.214064 515
## [5] {LUNCH BAG PINK POLKADOT} => {LUNCH BAG RED RETROSPOT} 0.02863083 0.5446328 7.770248 482
## [6] {LUNCH BAG RED RETROSPOT} => {LUNCH BAG PINK POLKADOT} 0.02863083 0.4084746 7.770248 482
## [7] {LUNCH BAG BLACK SKULL.} => {LUNCH BAG RED RETROSPOT} 0.02827443 0.4783920 6.825194 476
## [8] {LUNCH BAG RED RETROSPOT} => {LUNCH BAG BLACK SKULL.} 0.02827443 0.4033898 6.825194 476
## [9] {JUMBO BAG PINK POLKADOT} => {JUMBO BAG RED RETROSPOT} 0.02684883 0.5955204 7.641453 452
## [10] {PINK REGENCY TEACUP AND SAUCER} => {GREEN REGENCY TEACUP AND SAUCER} 0.02625483 0.8215613 20.643261 442Si ordenem per support, veiem les regles que més cops s’han produït a tot el conjunt. La regla comprar {GREEN REGENCY TEACUP AND SAUCER} => {ROSES REGENCY TEACUP AND SAUCER } s’ha produït més de 517 cops. Aquestes regles tenen a més una probabilitat condicionada (confidence) relativament alta > 40%, algunes amb més d’un 80%.
El que queda força clar que els productes són similars i d’aquí el patró de compra.
Ara executarem el mateix algoritme baixant el support i pujan el confidence.
rules2 <- apriori(transations, parameter = list(support = 0.005, confidence = 0.9))## Apriori
##
## Parameter specification:
## confidence minval smax arem aval originalSupport maxtime support minlen
## 0.9 0.1 1 none FALSE TRUE 5 0.005 1
## maxlen target ext
## 10 rules FALSE
##
## Algorithmic control:
## filter tree heap memopt load sort verbose
## 0.1 TRUE TRUE FALSE TRUE 2 TRUE
##
## Absolute minimum support count: 84
##
## set item appearances ...[0 item(s)] done [0.00s].
## set transactions ...[3578 item(s), 16835 transaction(s)] done [0.30s].
## sorting and recoding items ... [1161 item(s)] done [0.01s].
## creating transaction tree ... done [0.01s].
## checking subsets of size 1 2 3 4 5 6 done [0.07s].
## writing ... [221 rule(s)] done [0.00s].
## creating S4 object ... done [0.01s].inspect(head(sort(rules2, by = "confidence"), 10))## lhs rhs support confidence lift count
## [1] {HERB MARKER BASIL,
## HERB MARKER CHIVES ,
## HERB MARKER MINT,
## HERB MARKER ROSEMARY,
## HERB MARKER THYME} => {HERB MARKER PARSLEY} 0.007722008 0.9923664 90.79613 130
## [2] {GREEN REGENCY TEACUP AND SAUCER,
## REGENCY TEA PLATE PINK,
## REGENCY TEA PLATE ROSES } => {REGENCY TEA PLATE GREEN } 0.005524206 0.9893617 70.57587 93
## [3] {GREEN REGENCY TEACUP AND SAUCER,
## REGENCY TEA PLATE PINK,
## REGENCY TEA PLATE ROSES ,
## ROSES REGENCY TEACUP AND SAUCER } => {REGENCY TEA PLATE GREEN } 0.005227205 0.9887640 70.53323 88
## [4] {HERB MARKER CHIVES ,
## HERB MARKER MINT,
## HERB MARKER ROSEMARY,
## HERB MARKER THYME} => {HERB MARKER PARSLEY} 0.008197208 0.9857143 90.18750 138
## [5] {HERB MARKER BASIL,
## HERB MARKER CHIVES ,
## HERB MARKER ROSEMARY,
## HERB MARKER THYME} => {HERB MARKER PARSLEY} 0.007959608 0.9852941 90.14906 134
## [6] {HERB MARKER BASIL,
## HERB MARKER CHIVES ,
## HERB MARKER PARSLEY,
## HERB MARKER ROSEMARY} => {HERB MARKER THYME} 0.007959608 0.9852941 91.13971 134
## [7] {HERB MARKER BASIL,
## HERB MARKER CHIVES ,
## HERB MARKER MINT,
## HERB MARKER THYME} => {HERB MARKER PARSLEY} 0.007840808 0.9850746 90.12897 132
## [8] {HERB MARKER BASIL,
## HERB MARKER CHIVES ,
## HERB MARKER MINT,
## HERB MARKER ROSEMARY} => {HERB MARKER PARSLEY} 0.007840808 0.9850746 90.12897 132
## [9] {HERB MARKER BASIL,
## HERB MARKER CHIVES ,
## HERB MARKER MINT,
## HERB MARKER PARSLEY,
## HERB MARKER THYME} => {HERB MARKER ROSEMARY} 0.007722008 0.9848485 88.66270 130
## [10] {HERB MARKER BASIL,
## HERB MARKER CHIVES ,
## HERB MARKER MINT,
## HERB MARKER PARSLEY,
## HERB MARKER ROSEMARY} => {HERB MARKER THYME} 0.007722008 0.9848485 91.09848 130Ara hem modificat l’algoritme per a produir regles amb un confidence molt alt, encara que no es produeixin gaires cops en el total del dataset(support més baix). Tenim similaritat de productes, però encara així seria una bona forma de treure campanyes de marketing.
inspect(head(sort(rules2, by = "support"), 10))## lhs rhs support confidence lift count
## [1] {POPPY'S PLAYHOUSE BEDROOM ,
## POPPY'S PLAYHOUSE LIVINGROOM } => {POPPY'S PLAYHOUSE KITCHEN} 0.01069201 0.9000000 44.56324 180
## [2] {HERB MARKER THYME} => {HERB MARKER ROSEMARY} 0.01021681 0.9450549 85.08021 172
## [3] {HERB MARKER ROSEMARY} => {HERB MARKER THYME} 0.01021681 0.9197861 85.08021 172
## [4] {HERB MARKER PARSLEY} => {HERB MARKER ROSEMARY} 0.01015741 0.9293478 83.66615 171
## [5] {HERB MARKER ROSEMARY} => {HERB MARKER PARSLEY} 0.01015741 0.9144385 83.66615 171
## [6] {HERB MARKER MINT} => {HERB MARKER ROSEMARY} 0.01009801 0.9042553 81.40716 170
## [7] {HERB MARKER ROSEMARY} => {HERB MARKER MINT} 0.01009801 0.9090909 81.40716 170
## [8] {HERB MARKER BASIL} => {HERB MARKER ROSEMARY} 0.01003861 0.9184783 82.68760 169
## [9] {HERB MARKER ROSEMARY} => {HERB MARKER BASIL} 0.01003861 0.9037433 82.68760 169
## [10] {HERB MARKER PARSLEY} => {HERB MARKER MINT} 0.00997921 0.9130435 81.76110 168El mateix ordenat per support. Els lifts són altíssims, per sobre d’1 ens donaríem per satisfets.
inspect(head(sort(rules2, by = "lift"), 10))## lhs rhs support confidence lift count
## [1] {DOLLY GIRL CHILDRENS CUP,
## SPACEBOY CHILDRENS BOWL} => {DOLLY GIRL CHILDRENS BOWL} 0.005405405 0.9680851 120.72380 91
## [2] {DOLLY GIRL CHILDRENS BOWL,
## SPACEBOY CHILDRENS CUP} => {DOLLY GIRL CHILDRENS CUP} 0.005049005 0.9550562 116.50993 85
## [3] {DOLLY GIRL CHILDRENS BOWL,
## SPACEBOY CHILDRENS CUP} => {SPACEBOY CHILDRENS BOWL} 0.005049005 0.9550562 98.64031 85
## [4] {HERB MARKER BASIL,
## HERB MARKER MINT,
## HERB MARKER PARSLEY,
## HERB MARKER ROSEMARY,
## HERB MARKER THYME} => {HERB MARKER CHIVES } 0.007722008 0.9219858 95.81254 130
## [5] {HERB MARKER MINT,
## HERB MARKER PARSLEY,
## HERB MARKER ROSEMARY,
## HERB MARKER THYME} => {HERB MARKER CHIVES } 0.008197208 0.9139073 94.97302 138
## [6] {HERB MARKER BASIL,
## HERB MARKER MINT,
## HERB MARKER PARSLEY,
## HERB MARKER THYME} => {HERB MARKER CHIVES } 0.007840808 0.9103448 94.60281 132
## [7] {HERB MARKER BASIL,
## HERB MARKER PARSLEY,
## HERB MARKER ROSEMARY,
## HERB MARKER THYME} => {HERB MARKER CHIVES } 0.007959608 0.9054054 94.08951 134
## [8] {HERB MARKER BASIL,
## HERB MARKER MINT,
## HERB MARKER PARSLEY,
## HERB MARKER ROSEMARY} => {HERB MARKER CHIVES } 0.007840808 0.9041096 93.95485 132
## [9] {HERB MARKER MINT,
## HERB MARKER PARSLEY,
## HERB MARKER THYME} => {HERB MARKER CHIVES } 0.008375408 0.9038462 93.92747 141
## [10] {HERB MARKER BASIL,
## HERB MARKER CHIVES ,
## HERB MARKER PARSLEY,
## HERB MARKER ROSEMARY} => {HERB MARKER THYME} 0.007959608 0.9852941 91.13971 134La conclusió que podem treure d’aquest estudi és que a partir de l’algoritme d’associació podem veure patrons de conducta en les compres que es repeteixen i podem trobar associacions amb una fiabilitat molt alta.
Aquest algoritme primerament fixa un nombre de clusters (o centres de clusters) i a partir d’aquest nombre construeix els grups. Per a executar el mètode de k-means o centroides partim de la base que no coneixem el nombre òptim de clústers. Provarem amb diversos centres (del 2 al 10) per veure quina aproximació és millor.
Com funciona Kmeans?
Havíem preparat el model per fer dos estudis: de productes com de clients. Començarem pels clients.
customer.data <- read.csv('retail-customer.csv', sep = ',')
summary(customer.data)## NumInvoices NumCancellations TotalCost NumStocks
## Min. : 0.00 Min. : 0.0000 Min. : -4287.6 Min. : 0.00
## 1st Qu.: 1.00 1st Qu.: 0.0000 1st Qu.: 293.4 1st Qu.: 15.00
## Median : 2.00 Median : 0.0000 Median : 648.1 Median : 35.00
## Mean : 4.24 Mean : 0.8358 Mean : 1898.5 Mean : 61.03
## 3rd Qu.: 5.00 3rd Qu.: 1.0000 3rd Qu.: 1611.7 3rd Qu.: 77.00
## Max. :210.00 Max. :47.0000 Max. :279489.0 Max. :1787.00
## Days DiffInvoices
## Min. : 0.00 Min. : -6.000
## 1st Qu.: 16.10 1st Qu.: 1.000
## Median : 49.90 Median : 2.000
## Mean : 91.59 Mean : 3.404
## 3rd Qu.:142.93 3rd Qu.: 4.000
## Max. :373.10 Max. :196.000Anem a fer l’estudi amb 4 variables:
# filtrem els camps a treballar
customer.data.clean <- customer.data[,c('Days','NumInvoices','TotalCost','NumStocks')]
summary(customer.data.clean)## Days NumInvoices TotalCost NumStocks
## Min. : 0.00 Min. : 0.00 Min. : -4287.6 Min. : 0.00
## 1st Qu.: 16.10 1st Qu.: 1.00 1st Qu.: 293.4 1st Qu.: 15.00
## Median : 49.90 Median : 2.00 Median : 648.1 Median : 35.00
## Mean : 91.59 Mean : 4.24 Mean : 1898.5 Mean : 61.03
## 3rd Qu.:142.93 3rd Qu.: 5.00 3rd Qu.: 1611.7 3rd Qu.: 77.00
## Max. :373.10 Max. :210.00 Max. :279489.0 Max. :1787.00Eliminarem els outliers de totalCost, ja que crec que poden desvirtuar la realitat. Per a trobar valors extrems anem a aplicar la idea dels IQR (interquartile ranges):
Per a una determinada variable contínua, els outliers són aquelles observacions que es troben fora de 1.5 * IQR, on IQR, el “Inter Quartile Range” és la diferència entre el Q3 i el Q1:
# funció per treure els límits dels valors atípics:
outliersLimits <- function(x) {
limits <- c("above", "under")
limits$above <- quantile(x, 0.75, type=6) + 1.5 * (quantile(x, 0.75, type=6) - quantile(x, 0.25, type=6))
limits$under <- quantile(x, 0.25, type=6) - 1.5 * (quantile(x, 0.75, type=6) - quantile(x, 0.25, type=6))
return(limits)
}
limits <- outliersLimits(customer.data.clean$TotalCost)
# treiem aquest valors atípics de total cost. Per el límit superiot anem a agafar aquells que tenen una despesa superior a 0
customer.data.clean <-subset(customer.data.clean, TotalCost > 0)
customer.data.clean <-subset(customer.data.clean, TotalCost <= limits$above)
summary(customer.data.clean)## Days NumInvoices TotalCost NumStocks
## Min. : 0.00 Min. : 1.000 Min. : 2.9 Min. : 1.00
## 1st Qu.: 20.80 1st Qu.: 1.000 1st Qu.: 275.3 1st Qu.: 15.00
## Median : 56.10 Median : 2.000 Median : 581.8 Median : 31.00
## Mean : 97.27 Mean : 2.923 Mean : 869.3 Mean : 48.75
## 3rd Qu.:155.00 3rd Qu.: 4.000 3rd Qu.:1223.1 3rd Qu.: 65.00
## Max. :373.10 Max. :39.000 Max. :3580.1 Max. :580.00Com la magnitud dels valors difereix notablement entre les variables, les escalem abans d’aplicar el clustering.
df1 <- scale(customer.data.clean)Mostrem en una gràfica els valors de la silueta mitjana de cada prova per comprovar quin nombre de clusters és el millor. Com més alt és el valor de la silueta, millor
avg_sil <- function(k) {
km.res <- kmeans(df1, centers = k, nstart = 25, iter.max=25)
ss <- silhouette(km.res$cluster, dist(df1))
mean(ss[, 3])
}
k.values <- 2:12
avg_sil_values <- map_dbl(k.values, avg_sil)
plot(k.values, avg_sil_values,
type = "b", pch = 19, frame = FALSE,
xlab = "Nombre de clusters K",
ylab = "Avg. Silhouettes")
avg_sil_values## [1] 0.4549749 0.4022606 0.3755009 0.3620711 0.3456050 0.3495608 0.3122243
## [8] 0.2937851 0.2958048 0.2938006 0.2971479El millor valor obtingut és de 2 clusters o tipus de clients, però del 3 al 7 estan relativament molt a prop quant a valor de Avg. Silhouettes.
Podem executar el mateix procés amb la funció fviz_nbclust
fviz_nbclust(df1, kmeans, method = "silhouette")
Res de nou. Surten 2 clusters. I del 3 al 7, bons valors també.
Un altra forma d’avaluar quin és el millor nombre de clústers és considerar el millor model amb aquestes condicions:
Ofereix la menor suma dels quadrats de les distàncies dels punts de cada grup respecte al seu centre (withinss). Això vol dir que els elements dels grups estan junts.
La separació més gran entre centres de grups (betweenss). Els grups están separats.
És una idea conceptualment similar a la silueta. Una manera comú de fer la selecció del nombre de clústers consisteix a aplicar el mètode elbow (colze). Es fa una selecció del nombre de clústers a partir de la inspecció de la gràfica que s’obté l’iterar amb el mateix conjunt de dades per a distints valors del nombre de clústers. Es seleccionarà el valor que es troba en el colze de la corba.
wss <- function(k) {
kmeans(df1, k, nstart = 25, iter.max = 25 )$tot.withinss
}
# fem fins a 15 iteracions
k.values <- 1:15
wss_values <- map_dbl(k.values, wss)
plot(k.values, wss_values,
type="b", pch = 19, frame = FALSE,
xlab="Nombre de clusters K",
ylab="Total within-clusters sum of squares")
# la funció fviz_nbclust
fviz_nbclust(df1, kmeans, method = "wss")
A partir de la corba obtinguda podem veure com a mesura que s’augmenta la quantitat de centroides, el valor de “tot.tot.withinss” disminueix, ja que la separació entre elements dels clusters és inferior. La idea és trobar un “colze”. El colze es troba on ja no es produeixen variacions importants en augmentar ‘Nombre de clusters’. El valor que jo dedueixo és el 3, és on tenim el colze.
També es pot fer servir la funció kmeansruns del paquet fpc que executarà l’algoritme kmeans com un conjunt de valors i selecciona el valor del nombre de clústers que millor funcioni d’acord amb dos criteris: la silueta mitja (asw) i Calinski-Harabasz (“ch”).
fit_ch <- kmeansruns(df1, krange = 1:10, criterion = "ch", iter.max=25, runs=25)
plot(1:10,fit_ch$crit,type="o",col="blue",pch=0,xlab="Nombre de clústers",ylab="Criteri Calinski-Harabasz")
Amb aquest criteri estem obtenint 3 clusters també.
fit_asw <- kmeansruns(df1, krange = 1:10, criterion = "asw", iter.max=25, runs=25)
plot(1:10,fit_asw$crit,type="o",col="blue",pch=0,xlab="Nombre de clústers",ylab="Criteri silueta mitja")
Realitzem comparacions visuals amb el nombre de clusters que hem trobat en la majoria d’estudis (2 i 3 clusters)
# realitzem una divisió en 2 clusters
clusters <- kmeans(df1, centers = 2, iter.max = 25, nstart = 25)
clusters$centers## Days NumInvoices TotalCost NumStocks
## 1 -0.6280280 1.2333845 1.3431643 1.1268210
## 2 0.2151878 -0.4226074 -0.4602225 -0.3860945fviz_cluster(clusters, data = df1, geom = "point", repel = TRUE)
El mapa ens diferencia 2 grups de clients força clars.
fviz_cluster(clusters, data = df1, geom = "point", choose.vars = c("Days", "NumInvoices"), repel = TRUE, ggtheme = theme_minimal())
Si mostrem el nombre de factures x dies. Veiem una sèrie de clients potencialment molt bons, on han comprat relativament fa poc i tenen moltes factures. Són clients habituals.
fviz_cluster(clusters, data = df1, geom = "point", choose.vars = c("Days", "TotalCost"), repel = TRUE, ggtheme = theme_minimal())
En aquesta mostra tenim els millors clients en el grupuscle superior-esquerra. Més despesa i menys dies des de la darrera compra.
fviz_cluster(clusters, data = df1, geom = "point", choose.vars = c("NumInvoices", "TotalCost"), repel = TRUE, ggtheme = theme_minimal())
I en l’estudi de les dues variables, nombre de factures x cost. Diferenciem clients que gasten molt amb poques factures (part superior - esquerra) i altres que gasten força amb més factures (part superior - dreta).
Fem l’estudi visual amb 3 clusters:
# realitzem una divisió en 3 clusters
clusters <- kmeans(df1, centers = 3, iter.max = 25, nstart = 25)
clusters$centers## Days NumInvoices TotalCost NumStocks
## 1 -0.4628347 -0.2483648 -0.2928991 -0.2352042
## 2 1.5221179 -0.5952659 -0.6190990 -0.5331632
## 3 -0.6343581 1.4406513 1.5944061 1.3257665fviz_cluster(clusters, data = df1, geom = "point", repel = TRUE)
Un dels dos clusters que teníem s’ha dividit en dos (cluster 1 i 2).
set.seed(123)
fviz_cluster(clusters, data = df1, geom = "point", choose.vars = c("Days", "NumInvoices"), repel = TRUE, ggtheme = theme_minimal())
Mostrem el nombre de factures x dies:
set.seed(123)
fviz_cluster(clusters, data = df1, geom = "point", choose.vars = c("Days", "TotalCost"), repel = TRUE, ggtheme = theme_minimal())
En aquesta mostra tenim, Days vs Cost:
set.seed(123)
fviz_cluster(clusters, data = df1, geom = "point", choose.vars = c("NumInvoices", "TotalCost"), repel = TRUE, ggtheme = theme_minimal())
I en l’estudi de les dues variables, nombre de factures x cost:
En l’estudi hem optat per a realitzar les gràfiques visuals amb 2 i 3 clusters. De tota manera, els primers algoritmes ens estaven donant possibilitat de més agrupacions fins a un total de 7.
Ara farem un estudi similar amb els productes.
retail.product.data <- read.csv('retail-product.csv', sep = ',')
summary(retail.product.data)## UnitPrice TotalStockUnits Days TotalStockPrice
## Min. : 0.000 Min. :-1460.0 Min. : 0.00 Min. : -86541.1
## 1st Qu.: 1.010 1st Qu.: 67.0 1st Qu.: 0.90 1st Qu.: 142.7
## Median : 1.950 Median : 414.5 Median : 3.10 Median : 750.7
## Mean : 3.818 Mean : 1395.3 Mean : 46.06 Mean : 2861.7
## 3rd Qu.: 3.880 3rd Qu.: 1437.0 3rd Qu.: 32.02 3rd Qu.: 2373.4
## Max. :744.150 Max. :53215.0 Max. :373.10 Max. :1064825.1Anem a fer l’estudi amb 3 variables:
product.data.clean <- retail.product.data[,c('UnitPrice', 'Days','TotalStockUnits')]
summary(product.data.clean)## UnitPrice Days TotalStockUnits
## Min. : 0.000 Min. : 0.00 Min. :-1460.0
## 1st Qu.: 1.010 1st Qu.: 0.90 1st Qu.: 67.0
## Median : 1.950 Median : 3.10 Median : 414.5
## Mean : 3.818 Mean : 46.06 Mean : 1395.3
## 3rd Qu.: 3.880 3rd Qu.: 32.02 3rd Qu.: 1437.0
## Max. :744.150 Max. :373.10 Max. :53215.0# eliminem valors negatius o igual a 0
product.data.clean <- subset(product.data.clean, product.data.clean$TotalStockUnits > 0)
product.data.clean <- subset(product.data.clean, product.data.clean$UnitPrice > 0)
product.data.clean <- subset(product.data.clean, product.data.clean$Days > 0)
summary(product.data.clean)## UnitPrice Days TotalStockUnits
## Min. : 0.040 Min. : 0.10 Min. : 1
## 1st Qu.: 0.990 1st Qu.: 1.00 1st Qu.: 65
## Median : 1.950 Median : 3.80 Median : 371
## Mean : 3.699 Mean : 46.63 Mean : 1227
## 3rd Qu.: 3.840 3rd Qu.: 34.90 3rd Qu.: 1264
## Max. :744.150 Max. :373.10 Max. :53215Eliminarem els outliers de total cost, ja que crec que poden desvirtuar la realitat. Funció per treure els límits dels valors atípics:
outliersLimits <- function(x) {
limits <- c("above", "under")
limits$above <- quantile(x, 0.75, type=6) + 1.5 * (quantile(x, 0.75, type=6) - quantile(x, 0.25, type=6))
limits$under <- quantile(x, 0.25, type=6) - 1.5 * (quantile(x, 0.75, type=6) - quantile(x, 0.25, type=6))
return(limits)
}
limits <- outliersLimits(product.data.clean$UnitPrice)
# treiem aquest valors atípics de total cost
product.data.clean <-subset(product.data.clean, UnitPrice <= limits$above)
summary(product.data.clean)## UnitPrice Days TotalStockUnits
## Min. :0.040 Min. : 0.10 Min. : 1
## 1st Qu.:0.870 1st Qu.: 1.00 1st Qu.: 75
## Median :1.660 Median : 3.80 Median : 426
## Mean :2.327 Mean : 46.23 Mean : 1306
## 3rd Qu.:3.150 3rd Qu.: 32.90 3rd Qu.: 1375
## Max. :8.110 Max. :373.10 Max. :53215Com que no volem que l’algoritme de clúster no depengui d’una unitat variable arbitrària, comencem escalant / estandarditzant les dades mitjançant la funció scale.
df2 <- scale(product.data.clean)Anem a fer servir el mètode de la silueta per a trobar el millor nombre de clusters. En resum, l’enfocament mitjà de la silueta mesura la qualitat d’una agrupació. És a dir, determina el bé que cada objecte es troba dins del seu clúster. Una gran amplada mitjana de silueta indica una bona agrupació. El mètode de silueta mitjana calcula la silueta mitjana d’observacions per a diferents valors de k.
avg_sil <- function(k) {
km.res <- kmeans(df2, centers = k, nstart = 25, iter.max=25)
ss <- silhouette(km.res$cluster, dist(df2))
mean(ss[, 3])
}
k.values <- 2:15
avg_sil_values <- map_dbl(k.values, avg_sil)
plot(k.values, avg_sil_values,
type = "b", pch = 19, frame = FALSE,
xlab = "Nombre de clusters K",
ylab = "Avg. Silhouettes")
Segons el mètode de la silueta obtenim el millor nombre de clusters entre 4 i 5 tipus de productes diferents.
Podem executar el mateix procés amb la funció fviz_nbclust
fviz_nbclust(df2, kmeans, method = "silhouette")
4 és el nombre de clusters, seguit del 5 altre cop.
# function to compute total within-cluster sum of square
wss <- function(k) {
kmeans(df2, k, nstart = 25,iter.max=25)$tot.withinss
}
# Compute and plot wss for k = 1 to k = 15
k.values <- 1:15
# extract wss for 2-15 clusters
wss_values <- map_dbl(k.values, wss)
plot(k.values, wss_values,
type="b", pch = 19, frame = FALSE,
xlab="Number of clusters K",
ylab="Total within-clusters sum of squares")
Tenim el colze on ja no es produeixen variacions importants en augmentar ‘Nombre de clusters’. El valor que jo dedueixo és entre 4 i 5, altre cop.
Afortunadament, aquest procés per calcular el “mètode del colze” i la silueta s’ha agrupat en una única funció (fviz_nbclust):
fviz_nbclust(df2, kmeans, method = "wss")
No hi ha molt dubte. Tenim entre 4 tipus de productes diferents segons dies des de la darrera compra, unitats comprades i cost del producte.
Anem a fer servir la funció kmeansruns del paquet fpc que executarà l’algoritme kmeans com un conjunt de valors i selecciona el valor del nombre de clústers que millor funcioni d’acord amb dos criteris: la silueta mitja (asw) i Calinski-Harabasz (“ch”).
fit_asw <- kmeansruns(df2, krange = 2:20, criterion = "asw", iter.max=25, runs=10)
plot(1:20,fit_asw$crit,type="o",col="blue",pch=0,xlab="Nombre de clústers",ylab="Criteri silueta mitja")
Tornem a tenir 4 o 5 clusters. Farem estudis visuals amb aquests valors, ja que els resultats són molt més clars que amb els customers.
# realitzem una divisió en 4 clusters
clusters4 <- kmeans(df2, centers = 4, iter.max=25, nstart = 25)
clusters4$centers## UnitPrice Days TotalStockUnits
## 1 -0.56788319 -0.5063242 3.62724749
## 2 1.51580309 -0.2952085 -0.27705397
## 3 -0.47872208 -0.3864604 -0.05177185
## 4 -0.03026779 2.2455655 -0.42090266fviz_cluster(clusters4, data = df2, geom = "point", repel = TRUE)
set.seed(123)
fviz_cluster(clusters4, data = df2, geom = "point", choose.vars = c("UnitPrice", "TotalStockUnits"), repel = TRUE, ggtheme = theme_minimal())
Unitats venudes totals x Preu unitari:
set.seed(123)
fviz_cluster(clusters4, data = df2, geom = "point", choose.vars = c("UnitPrice", "Days"), repel = TRUE, ggtheme = theme_minimal())
Preu unitari x dies des de la darrera compra:
set.seed(123)
fviz_cluster(clusters4, data = df2, geom = "point", choose.vars = c("TotalStockUnits", "Days"), repel = TRUE, ggtheme = theme_minimal())
Dies des de la darrera compra x Unitats venudes:
# realitzem una divisió en 5 clusters
clusters5 <- kmeans(df2, centers = 5, iter.max=25, nstart = 25)
clusters5$centers## UnitPrice Days TotalStockUnits
## 1 -0.50342869 -0.5161059 9.7281381
## 2 -0.47372940 -0.3825041 -0.1214569
## 3 1.51933002 -0.2939350 -0.2847242
## 4 -0.55164654 -0.4897042 2.2492208
## 5 -0.03026779 2.2455655 -0.4209027summary(df2)## UnitPrice Days TotalStockUnits
## Min. :-1.2337 Min. :-0.5265 Min. :-0.46919
## 1st Qu.:-0.7859 1st Qu.:-0.5163 1st Qu.:-0.44259
## Median :-0.3597 Median :-0.4843 Median :-0.31645
## Mean : 0.0000 Mean : 0.0000 Mean : 0.00000
## 3rd Qu.: 0.4442 3rd Qu.:-0.1521 3rd Qu.: 0.02471
## Max. : 3.1201 Max. : 3.7312 Max. :18.65567fviz_cluster(clusters5, data = df2, geom = "point", repel = TRUE)
Realment no ens aporta gaire generar un altre cluster. Amb 4 clusters ens quedem.
Funciona amb l’algoritme dels veïns més propers (k-nearest neighbours) de la següent forma:
Assignem cada element al seu propi clúster, de manera que si té N elements, ara té N clústers, cadascun amb un sol element. Deixem que les distàncies (similitud) entre els clusters siguin iguals a les distàncies (similitud) entre els elements que contenen.
Es troba el parell més proper (més similar) de clústers i es combinen en un únic clúster, de manera que ara tenim un clúster menys.
Calculem distàncies (similitud) entre el nou clúster i cadascun dels clústers antics.
Es repeteix els passos 2 i 3 fins que tots els elements siguin agrupats en un sol grup de mida N.
De formes de mesurar les distàncies entre grups tenim de diversos tipus. Provarem només: “average”, “complete” i “ward”.
Fem servir el dataframe: ‘customer.data.clean’ que hem generat a Clustering amb K-means o centroides:
summary(customer.data.clean)## Days NumInvoices TotalCost NumStocks
## Min. : 0.00 Min. : 1.000 Min. : 2.9 Min. : 1.00
## 1st Qu.: 20.80 1st Qu.: 1.000 1st Qu.: 275.3 1st Qu.: 15.00
## Median : 56.10 Median : 2.000 Median : 581.8 Median : 31.00
## Mean : 97.27 Mean : 2.923 Mean : 869.3 Mean : 48.75
## 3rd Qu.:155.00 3rd Qu.: 4.000 3rd Qu.:1223.1 3rd Qu.: 65.00
## Max. :373.10 Max. :39.000 Max. :3580.1 Max. :580.00Abans de tot, normalitzem les dades i calculem la distància euclidiana entre els elements.
## creem la funció per normalitzar
ni <-function (x) {(x -min (x)) / (max (x) -min (x))}
customer.data.norm <- as.data.frame(lapply(customer.data.clean[,1:4], ni))
# primer fem la matriu de la distància euclideana dels elements de la taula de llavors.
customer.data.euc <- dist(customer.data.norm, method = "euclidean")Estratègia de la distància màxima o similitud mínima. En aquest mètode, també conegut (complete linkage), es considera que la distància o similitud entre dos clústers cal mesurar-la atenent a elements més dispars, és a dir, la distància o similitud entre clústers ve dada, respectivament, per la màxima distància (o mínima similitud) entre components dels clústers.
# fem servir hclust per a realitzar les agrupacions
hc1 <- hclust(customer.data.euc, method = "complete")
plot(hc1,
main = "Complete linkage",
xlab = "Customers",
ylab = "",
cex = 0.04,
sub = "")
abline(h = 0.78, col = "red")
abline(h = 1.1, col = "blue")
El dendograma ens està mostrant les distàncies euclidianes entre els diferents elements (customers) de la mostra i la separació gràfica conformes en van agrupant els clusters. Les interpretació és subjectiva però agafariem 4 i 7 clusters.
En aquesta estratègia la distància, o similitud, de l’clúster Ci amb el Cj s’obté com la mitjana aritmètica entre la distància, o similitud, dels elements d’aquests clústers.
hc2 <- hclust(customer.data.euc, method = "average")
plot(hc2,
main = "Average linkage",
xlab = "Customers",
ylab = "",
cex = 0.1,
sub = "")
abline(h = 0.4, col = "red")
abline(h = 0.55, col = "blue")
Similar, entre 3 i 7 grups de customers.
El mètode de Ward apunta a minimitzar la variància total dins de el grup. A cada pas, es fusionen el parell de clústers amb una distància mínima entre els clústers. En altres paraules, forma grups d’una manera que minimitza la pèrdua associada amb cada grup. Té molt bona efectivitat a l’hora d’agrupar.
hc3 <- hclust(customer.data.euc, method = "ward.D")
plot(hc3,
main = "Ward Method",
xlab = "Customers",
ylab = "",
cex = 0.1,
sub = "")
abline(h = 100, col = "red")
abline(h = 55, col = "blue")
Amb el mètode ward veiem més clarament 4 o 5 grups de clients. Els altres clusters es troben molt més allunyats en l’agrupació. Ens quedariem amb 4 o 5 grups de customers.
Anem a reduir la quantitat d’elements de la mostra i apliquem algoritme ward
set.seed(123)
sample <- customer.data.norm %>% sample_frac(0.10)
sample.euc <- dist(sample, method = "euclidean")
hc3 <- hclust(sample.euc, method = "ward.D")
plot(hc3,
main = "Ward Method",
xlab = "Products",
ylab = "",
cex = 0.1,
sub = "")
abline(h = 12, col = "red")
abline(h = 7, col = "green")
Amb menys mostra estem treient entre 3 i 5 clusters.
Dels 3 mètodes que hem fet servir el que dona una més bona interpretació és ward on podem veure entre 3 i 5 clusters.
Fem servir el dataframe: ‘product.data.clean’ que hem generat a Clustering amb K-means o centroides:
summary(product.data.clean)## UnitPrice Days TotalStockUnits
## Min. :0.040 Min. : 0.10 Min. : 1
## 1st Qu.:0.870 1st Qu.: 1.00 1st Qu.: 75
## Median :1.660 Median : 3.80 Median : 426
## Mean :2.327 Mean : 46.23 Mean : 1306
## 3rd Qu.:3.150 3rd Qu.: 32.90 3rd Qu.: 1375
## Max. :8.110 Max. :373.10 Max. :53215Abans de tot, normalitzem les dades i calculem la distància euclidiana entre els elements.
## creem la funció per normalitzar
ni <-function (x) {(x -min (x)) / (max (x) -min (x))}
product.data.norm <- as.data.frame(lapply(product.data.clean[,1:3], ni))
# primer fem la matriu de la distància euclideana dels elements de la taula de llavors.
product.data.euc <- dist(product.data.norm, method = "euclidean")Executarem directament amb el mètode ward on veiem que tenim 4 clusters diferenciats.
hc3 <- hclust(product.data.euc, method = "ward.D")
plot(hc3,
main = "Ward Method",
xlab = "Products",
ylab = "",
cex = 0.1,
sub = "")
abline(h = 70, col = "red")
abline(h = 50, col = "blue")
El dendograma ens està mostrant les distàncies euclidianes entre els diferents elements (productes) de la mostra. Mostra la separació dels elements i dels clusters que es van formant. Sembla distingir entre 4 i 5 clusters. La distància fins a 3 clusters ja és molt llunyana.
Reduïm la dimensionalitat de la mostra i apliquem algoritme ward.
set.seed(123)
sample <- product.data.norm %>% sample_frac(0.10)
sample.euc <- dist(sample, method = "euclidean")
hc3 <- hclust(sample.euc, method = "ward.D")
plot(hc3,
main = "Ward Method",
xlab = "Products",
ylab = "",
cex = 0.1,
sub = "")
abline(h = 8, col = "red")
abline(h = 6, col = "green")
Podríem fer entre 4 i 5 clusters. Resultat molt similar al K-means o centroides.
Els enfocaments dels mètodes probabilístics assumeixen una varietat de models de dades i apliquen l’estimació de màxima versemblança i els criteris de Bayes per identificar el model més probable i el nombre de grups.
Agafem el model name VVV, (multivariate mixture, ellipsoidal, varying volume, shape, and orientation). És el nostre cas.
sample <- scale(customer.data.clean, center=TRUE, scale=TRUE)
mclust <- Mclust(sample[,1:4], modelName = "VVV")
summary(mclust)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 3
## components:
##
## log-likelihood n df BIC ICL
## -13584.47 3895 44 -27532.71 -28305.64
##
## Clustering table:
## 1 2 3
## 1340 1044 1511plot(mclust, what = 'BIC', legendArgs = list(x = "bottomright", ncol = 10))
plot(mclust, what = "uncertainty")
Ens ha generat 3 grups.
Anem a veure amb variables 2 a 2: ### Days x TotalCost
mclust <- Mclust(sample[,c("Days","TotalCost")], modelName = "VVV")
summary(mclust)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 9
## components:
##
## log-likelihood n df BIC ICL
## -7188.852 3895 53 -14815.88 -17161.94
##
## Clustering table:
## 1 2 3 4 5 6 7 8 9
## 625 176 411 414 687 444 98 582 458plot(mclust, what = 'BIC', legendArgs = list(x = "bottomright", ncol = 5), modelName = "VVV")
plot(mclust, what = "uncertainty")
Amb la l’estudi de Days i total cost ens està mostrant fins a 9 grups diferents. De 4 a 9 grups el BIC (Bayesian Information Criterion) és molt similar.
mclust <- Mclust(sample[,c("Days","NumStocks")], modelName = "VVV")
summary(mclust)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 9
## components:
##
## log-likelihood n df BIC ICL
## -7123.542 3895 53 -14685.26 -17379.11
##
## Clustering table:
## 1 2 3 4 5 6 7 8 9
## 415 422 602 743 510 498 338 99 268plot(mclust, what = 'BIC', legendArgs = list(x = "bottomright", ncol = 5), modelName = "VVV")
plot(mclust, what = "uncertainty")
Days x NumStocks ens passa exactament el mateix.
Anem a fer un darrer estudi del BIC
BIC <- mclustBIC(sample)
plot(BIC)
summary(BIC)## Best BIC values:
## VEV,9 VEV,8 VEV,7
## BIC -20733.6 -21255.9952 -23332.892
## BIC diff 0.0 -522.3999 -2599.297mod1 <- Mclust(sample, x = BIC)
summary(mod1, parameters = TRUE)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VEV (ellipsoidal, equal shape) model with 9 components:
##
## log-likelihood n df BIC ICL
## -9912.088 3895 110 -20733.6 -21589.19
##
## Clustering table:
## 1 2 3 4 5 6 7 8 9
## 863 162 329 675 152 298 710 592 114
##
## Mixing probabilities:
## 1 2 3 4 5 6 7
## 0.21259164 0.04382844 0.08530769 0.16654909 0.04180480 0.08144189 0.18138048
## 8 9
## 0.15040885 0.03668712
##
## Means:
## [,1] [,2] [,3] [,4] [,5]
## Days -0.8216975 -0.09640612 0.003541239 -0.1265173 0.326200610
## NumInvoices 0.8623392 0.32540992 0.105643460 -0.3588188 0.005405116
## TotalCost 0.8399817 0.52682303 -0.069819382 -0.4124888 -0.443884545
## NumStocks 0.6317835 0.03190393 -0.119282024 -0.3745369 0.440404839
## [,6] [,7] [,8] [,9]
## Days -0.4134963 -0.3295983 1.7172975 0.5780150
## NumInvoices 1.1834385 -0.7475808 -0.7475808 0.1251766
## TotalCost 1.3064223 -0.7232285 -0.7645318 0.8538076
## NumStocks 1.2386544 -0.5842337 -0.6478521 0.5714852
##
## Variances:
## [,,1]
## Days NumInvoices TotalCost NumStocks
## Days 0.006662607 -0.04057397 -0.03581171 -0.02021258
## NumInvoices -0.040573972 1.68784994 1.26853213 0.42003817
## TotalCost -0.035811706 1.26853213 1.71647894 0.60608890
## NumStocks -0.020212585 0.42003817 0.60608890 1.34835971
## [,,2]
## Days NumInvoices TotalCost NumStocks
## Days 0.8626668 -0.2272808 -0.1330553 -0.2681446
## NumInvoices -0.2272808 0.5775007 0.4543723 0.4742551
## TotalCost -0.1330553 0.4543723 1.1673724 1.0738690
## NumStocks -0.2681446 0.4742551 1.0738690 1.0224531
## [,,3]
## Days NumInvoices TotalCost NumStocks
## Days 1.3401967 -0.25010132 -0.2553406 -0.24003980
## NumInvoices -0.2501013 0.24985736 0.0811263 0.05356166
## TotalCost -0.2553406 0.08112630 0.3645501 0.33102521
## NumStocks -0.2400398 0.05356166 0.3310252 0.30788137
## [,,4]
## Days NumInvoices TotalCost NumStocks
## Days 4.125584e-01 4.437192e-20 -1.033005e-02 -1.940712e-02
## NumInvoices 4.437192e-20 7.125119e-04 2.645351e-17 -1.673268e-18
## TotalCost -1.033005e-02 2.645351e-17 1.017483e-01 4.379301e-02
## NumStocks -1.940712e-02 -1.673268e-18 4.379301e-02 9.161723e-02
## [,,5]
## Days NumInvoices TotalCost NumStocks
## Days 2.7065908 -0.3722718 -0.3036071 -0.5224246
## NumInvoices -0.3722718 0.9301467 0.4021086 0.1401498
## TotalCost -0.3036071 0.4021086 0.2359166 0.2183954
## NumStocks -0.5224246 0.1401498 0.2183954 0.4710598
## [,,6]
## Days NumInvoices TotalCost NumStocks
## Days 0.1367154 -0.6524757 -0.1597707 -0.9568664
## NumInvoices -0.6524757 7.3456921 2.0439347 3.8724098
## TotalCost -0.1597707 2.0439347 2.4725438 1.9786256
## NumStocks -0.9568664 3.8724098 1.9786256 9.4840355
## [,,7]
## Days NumInvoices TotalCost NumStocks
## Days 1.579166e-01 -2.116084e-19 1.602121e-03 -7.442771e-03
## NumInvoices -2.116084e-19 2.724293e-04 1.507792e-17 1.182417e-17
## TotalCost 1.602121e-03 1.507792e-17 3.116479e-02 1.558987e-02
## NumStocks -7.442771e-03 1.182417e-17 1.558987e-02 4.259392e-02
## [,,8]
## Days NumInvoices TotalCost NumStocks
## Days 1.738108e-01 -1.054055e-19 -1.566002e-03 -2.286618e-03
## NumInvoices -1.054055e-19 2.990525e-04 -7.976302e-17 4.282695e-17
## TotalCost -1.566002e-03 -7.976302e-17 3.644557e-02 1.784302e-02
## NumStocks -2.286618e-03 4.282695e-17 1.784302e-02 4.405941e-02
## [,,9]
## Days NumInvoices TotalCost NumStocks
## Days 2.429545 -1.2358349 -1.7067809 -1.634157
## NumInvoices -1.235835 0.6727232 0.8038324 1.139690
## TotalCost -1.706781 0.8038324 2.7577851 1.726544
## NumStocks -1.634157 1.1396904 1.7265439 5.955426sample <- scale(product.data.clean, center=TRUE, scale=TRUE)
mclust <- Mclust(sample[,1:3], modelName = "VVV")
summary(mclust)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 9
## components:
##
## log-likelihood n df BIC ICL
## -2062.095 3380 89 -4847.372 -6161.645
##
## Clustering table:
## 1 2 3 4 5 6 7 8 9
## 585 350 608 450 228 173 584 167 235plot(mclust, what = 'BIC', legendArgs = list(x = "bottomright", ncol = 5), modelName = "VVV")
plot(mclust, what = "uncertainty")
9 clusters en total. I de 6 a 9 clusters valors molt similars.
mclust <- Mclust(sample[,c("Days","UnitPrice")], modelName = "VVV")
summary(mclust)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 8
## components:
##
## log-likelihood n df BIC ICL
## -2072.217 3380 47 -4526.34 -5710.032
##
## Clustering table:
## 1 2 3 4 5 6 7 8
## 725 516 525 261 409 92 454 398plot(mclust, what = 'BIC', legendArgs = list(x = "bottomright", ncol = 5), modelName = "VVV")
plot(mclust, what = "uncertainty")
UnitPrice x Days ens dóna 9 clusters, pero de 4 a 8 és molt similar.
mclust <- Mclust(sample[,c("UnitPrice","TotalStockUnits")], modelName = "VVV")
summary(mclust)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 6
## components:
##
## log-likelihood n df BIC ICL
## -3923.13 3380 35 -8130.657 -9556.23
##
## Clustering table:
## 1 2 3 4 5 6
## 950 190 524 463 536 717plot(mclust, what = 'BIC', legendArgs = list(x = "bottomright", ncol = 5), modelName = "VVV")
plot(mclust, what = "uncertainty")
UnitPrice x TotalStockUnits de 7 a 9 clusters
mclust <- Mclust(sample[,c("TotalStockUnits","Days")], modelName = "VVV")
summary(mclust)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VVV (ellipsoidal, varying volume, shape, and orientation) model with 8
## components:
##
## log-likelihood n df BIC ICL
## 2519.101 3380 47 4656.297 3341.915
##
## Clustering table:
## 1 2 3 4 5 6 7 8
## 234 632 128 375 412 894 422 283plot(mclust, what = 'BIC', legendArgs = list(x = "bottomright", ncol = 5), modelName = "VVV")
plot(mclust, what = "uncertainty")
Anem a fer un altre estudi del BIC
BIC <- mclustBIC(sample)
plot(BIC)
summary(BIC)## Best BIC values:
## VVE,9 VVV,9 VVI,9
## BIC -3919.759 -4396.0839 -4483.5522
## BIC diff 0.000 -476.3252 -563.7936mod1 <- Mclust(sample, x = BIC)
summary(mod1, parameters = TRUE)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VVE (ellipsoidal, equal orientation) model with 9 components:
##
## log-likelihood n df BIC ICL
## -1695.796 3380 65 -3919.759 -5298.618
##
## Clustering table:
## 1 2 3 4 5 6 7 8 9
## 566 486 70 459 435 278 266 466 354
##
## Mixing probabilities:
## 1 2 3 4 5 6 7
## 0.14716265 0.14459761 0.02238242 0.14674063 0.12349818 0.09292232 0.08689564
## 8 9
## 0.12104182 0.11475873
##
## Means:
## [,1] [,2] [,3] [,4] [,5]
## UnitPrice -0.7315634 0.5734528 -0.7112881 0.8292693 0.2639911
## Days -0.5036560 1.3779237 0.1362377 -0.5063829 -0.4223195
## TotalStockUnits 0.1180646 -0.4652270 2.2871462 -0.2564182 -0.4207641
## [,6] [,7] [,8] [,9]
## UnitPrice -0.4865466 -0.2410860 -0.28897304 -0.1088588
## Days -0.3263965 -0.5200799 -0.51688709 1.1883656
## TotalStockUnits -0.1603981 1.7290312 -0.02533168 -0.3832369
##
## Variances:
## [,,1]
## UnitPrice Days TotalStockUnits
## UnitPrice 7.293574e-02 -5.268068e-05 8.367072e-05
## Days -5.268068e-05 2.517675e-04 6.261126e-05
## TotalStockUnits 8.367072e-05 6.261126e-05 1.812563e-01
## [,,2]
## UnitPrice Days TotalStockUnits
## UnitPrice 1.3439005606 0.0003799178 -1.037973e-03
## Days 0.0003799178 1.8682837538 -6.461269e-04
## TotalStockUnits -0.0010379728 -0.0006461269 1.305055e-05
## [,,3]
## UnitPrice Days TotalStockUnits
## UnitPrice 0.1581333816 0.0002837577 0.011152469
## Days 0.0002837577 0.5440980293 0.004858508
## TotalStockUnits 0.0111524691 0.0048585084 14.599448818
## [,,4]
## UnitPrice Days TotalStockUnits
## UnitPrice 1.2136916781 -8.803005e-04 -9.163063e-04
## Days -0.0008803005 2.385458e-04 9.856040e-06
## TotalStockUnits -0.0009163063 9.856040e-06 2.678426e-02
## [,,5]
## UnitPrice Days TotalStockUnits
## UnitPrice 0.8506246414 -6.132157e-04 -6.557828e-04
## Days -0.0006132157 5.343750e-03 -9.640934e-07
## TotalStockUnits -0.0006557828 -9.640934e-07 1.185384e-03
## [,,6]
## UnitPrice Days TotalStockUnits
## UnitPrice 0.2681562487 -1.782899e-04 -1.661928e-04
## Days -0.0001782899 2.238269e-02 1.067855e-05
## TotalStockUnits -0.0001661928 1.067855e-05 5.287550e-02
## [,,7]
## UnitPrice Days TotalStockUnits
## UnitPrice 0.6992368420 -5.067393e-04 0.0009117609
## Days -0.0005067393 3.346052e-05 0.0006501426
## TotalStockUnits 0.0009117609 6.501426e-04 1.8796373845
## [,,8]
## UnitPrice Days TotalStockUnits
## UnitPrice 0.3194221003 -2.316995e-04 -1.696151e-04
## Days -0.0002316995 1.040126e-06 3.463848e-05
## TotalStockUnits -0.0001696151 3.463848e-05 9.968597e-02
## [,,9]
## UnitPrice Days TotalStockUnits
## UnitPrice 0.8546688172 0.0001050217 -0.0006549083
## Days 0.0001050217 0.9998011482 -0.0003433872
## TotalStockUnits -0.0006549083 -0.0003433872 0.0066835656Els millors BIC (Bayesian Information Criterion) que tenim és de 9 clusters.
Partim de la base que les dades que tenim són de comportament de compra dels clients i de tipologia de productes:
Fem un resum dels resultats obtinguts amb els 3 tipus de fer clustering pels clients.
Hem decidit que teníem entre 2 i 3 clusters, però en algun dels estudis hem tingut fins a 7 possibles clusters.
Són els que ens han donat uns resultats més clars. Encara que amb enllaç complet i promig hem tingut de 3 a 7 clusters amb mètode ward hem pogut fixar millor el nombre de clusters a 4 o 5.
En aquest cas hem tingut més disparitat de dades. En el primer estudi ens ha donat 3 clusters i després 5,6 i 9. No m’ha generat molta fiabilitat l’estudi fet amb mètodes probabilistics.
De totes maneres hem vist en tots els estudis que hi havia un ventall una mica gran de possibles tipologies de clients, però si hagués de fer un estudi, seria entre 3 i 5 clusters.
Fem el mateix resum però pels productes.
En aquest cas hem tingut molt clarament una agrupació entre 4 i 5 clusters i més concretament en 4.
Mateixos resultats que amb K-Means, 4 i 5 clusters.
Un altre cop tenim valors molt dispars comparats amb els mètodes basats en les distàncies HAC i K-Means.
Per a crear el model de l’arbre de decisió anem a fer servir la llibreria C50. Aquesta llibreria implementa funcions per a generar arbres amb l’algoritme C5.0 alhora basat en l’algoritme ID3 de Quinlan. Com a millores, permet treballar amb variables categòriques i discretes i realitza la poda de forma automàtica.
El mètode ID3 tracta de trobar una partició que asseguri la màxima capacitat predictiva i la màxima homogeneïtat de les classes. Volem pocs atributs per a predir el màxim possible. Per a fer-ho, a l’escollir un atribut, calcula el guany obtingut, i el càlcul el fa amb la diferència de la informació d’una partició respecte al desordre (o entropia) que es genera a partir d’aquella elecció.
La millora en la poda la realitza calculant la taxa d’error de les prediccions de les fulles. Es realitza una estimació pessimista de la predicció (amb calcul de la distribucuó binomial), i si els nodes fills tenen una estimació pitjor que la dels pares s’eliminen. És postpoda, és a dir, es realitza la poda després de contruir l’arbre.
En aquest primer estudi anem a executar l’algoritme de predicció sense excloure cap dels atributs del dataset.
Prepararem les dades per a avaluar l’arbre de decisió. Dividim en conjunt de prova i conjunt d’entrenament. Un el fem servir per a construir el model de l’arbre de decisió i el de prova per avaluar quina precisió obtenim d’aquest arbre. Fem servir la proporció 0.7 entrenament, 0.3 d’avaluació. La variable a avaluar en aquest cas és la columna ‘y’.
# separem el 70% de conjunt d'entrenament del 30% per a avaluar
## obtenim primer una mostra del total de forma aleatoria (hem guardat físicament aquest random per poder reproduïr el mateix i el carreguem de disc)
# random.data <- sample(1:nrow(bank.data), 0.70 * nrow(bank.data))
# write.csv(random.data,"random-data.csv",row.names=FALSE)
bank.data <- read.csv('bank-clean-discret.csv', sep = ',')
summary(bank.data)## age job marital education
## Min. :1.000 blue-collar:8507 divorced: 4494 primary : 5926
## 1st Qu.:2.000 management :7896 married :22762 secondary:20981
## Median :3.000 technician :6615 single :10997 tertiary :11346
## Mean :3.046 admin. :4565
## 3rd Qu.:4.000 services :3708
## Max. :5.000 retired :1492
## (Other) :5470
## default balance housing loan contact
## no :37485 Min. :-1884.0 no :16220 no :31548 cellular :24919
## yes: 768 1st Qu.: 43.0 yes:22033 yes: 6705 telephone: 2169
## Median : 338.0 unknown :11165
## Mean : 623.6
## 3rd Qu.: 957.0
## Max. : 3388.0
##
## day month duration campaign
## Min. : 1.00 may :12161 Min. :1.000 Min. : 1.000
## 1st Qu.: 8.00 jul : 6141 1st Qu.:2.000 1st Qu.: 1.000
## Median :16.00 aug : 5320 Median :3.000 Median : 2.000
## Mean :15.78 jun : 4333 Mean :3.003 Mean : 2.776
## 3rd Qu.:21.00 nov : 2940 3rd Qu.:4.000 3rd Qu.: 3.000
## Max. :31.00 apr : 2450 Max. :5.000 Max. :58.000
## (Other): 4908
## pdays previous poutcome y
## Min. : -1.0 Min. : 0.0000 failure: 4099 no :34124
## 1st Qu.: -1.0 1st Qu.: 0.0000 other : 1551 yes: 4129
## Median : -1.0 Median : 0.0000 success: 1156
## Mean : 40.3 Mean : 0.5659 unknown:31447
## 3rd Qu.: -1.0 3rd Qu.: 0.0000
## Max. :871.0 Max. :275.0000
## random.data <- read.csv('random-data.csv', sep = ',')
random.data <- unlist(random.data, use.names=FALSE)
data.train <- bank.data[random.data,]
data.test <- bank.data[-random.data,]
# en les dades d'entrenament, separem les dades per a generar l'arbre de la dada objectiu (target)
train.x <- data.train[,1:16]
train.y <- data.train[,17]
# en les dades de test, separem les dades per a generar l'arbre de la dada objectiu (target)
test.x <- data.test[,1:16]
test.y <- data.test[,17]model1 <- C50::C5.0(train.x, train.y, rules=TRUE)
predicted.model1 <- predict(model1, test.x, type="class")
cat(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model1 == test.y) / length(predicted.model1)))## La precissió de l'arbre és del: 91.20 %summary(model1)##
## Call:
## C5.0.default(x = train.x, y = train.y, rules = TRUE)
##
##
## C5.0 [Release 2.07 GPL Edition] Tue Jun 16 22:30:08 2020
## -------------------------------
##
## Class specified by attribute `outcome'
##
## Read 26777 cases (17 attributes) from undefined.data
##
## Rules:
##
## Rule 1: (3852/22, lift 1.1)
## contact in {telephone, unknown}
## duration <= 2
## -> class no [0.994]
##
## Rule 2: (152, lift 1.1)
## housing = yes
## month = jan
## pdays > 155
## -> class no [0.994]
##
## Rule 3: (1168/10, lift 1.1)
## housing = no
## contact = cellular
## day > 7
## month = aug
## duration <= 2
## -> class no [0.991]
##
## Rule 4: (7043/79, lift 1.1)
## month in {apr, dec, jan, jul, jun, may}
## duration <= 2
## -> class no [0.989]
##
## Rule 5: (782/11, lift 1.1)
## contact = cellular
## month = nov
## duration <= 2
## -> class no [0.985]
##
## Rule 6: (9851/159, lift 1.1)
## duration <= 2
## previous <= 2
## poutcome in {failure, other, unknown}
## -> class no [0.984]
##
## Rule 7: (9138/155, lift 1.1)
## housing = yes
## duration <= 3
## -> class no [0.983]
##
## Rule 8: (4770/112, lift 1.1)
## day > 20
## duration <= 3
## pdays <= 385
## poutcome in {failure, other, unknown}
## -> class no [0.976]
##
## Rule 9: (14615/367, lift 1.1)
## balance <= 2174
## duration <= 3
## poutcome in {failure, other, unknown}
## -> class no [0.975]
##
## Rule 10: (7676/291, lift 1.1)
## contact = unknown
## month in {jan, jul, jun, may}
## poutcome in {failure, other, unknown}
## -> class no [0.962]
##
## Rule 11: (930/45, lift 1.1)
## housing = yes
## day > 17
## month = nov
## -> class no [0.951]
##
## Rule 12: (2896/141, lift 1.1)
## housing = yes
## campaign > 3
## previous <= 5
## -> class no [0.951]
##
## Rule 13: (356/20, lift 1.1)
## housing = yes
## day <= 4
## month = feb
## previous <= 3
## -> class no [0.941]
##
## Rule 14: (3947/241, lift 1.1)
## education = primary
## pdays <= 385
## previous <= 4
## poutcome in {failure, other, unknown}
## -> class no [0.939]
##
## Rule 15: (15016/970, lift 1.0)
## housing = yes
## pdays <= 385
## poutcome in {failure, other, unknown}
## -> class no [0.935]
##
## Rule 16: (18970/1323, lift 1.0)
## month in {aug, jan, jul, may, nov}
## pdays <= 385
## poutcome in {failure, other, unknown}
## -> class no [0.930]
##
## Rule 17: (1062/73, lift 1.0)
## marital in {divorced, married}
## day <= 9
## month in {apr, feb}
## pdays <= 385
## poutcome in {failure, other, unknown}
## -> class no [0.930]
##
## Rule 18: (708/49, lift 1.0)
## day <= 20
## month in {apr, feb}
## pdays <= 385
## poutcome = failure
## -> class no [0.930]
##
## Rule 19: (1357/99, lift 1.0)
## day <= 9
## month = feb
## pdays <= 385
## poutcome in {failure, other, unknown}
## -> class no [0.926]
##
## Rule 20: (15, lift 8.7)
## housing = no
## month = feb
## pdays <= 92
## poutcome = success
## -> class yes [0.941]
##
## Rule 21: (24/2, lift 8.1)
## age > 3
## housing = no
## contact = cellular
## month = feb
## poutcome = success
## -> class yes [0.885]
##
## Rule 22: (23/2, lift 8.1)
## marital = divorced
## day > 17
## duration > 2
## campaign <= 3
## poutcome = success
## -> class yes [0.880]
##
## Rule 23: (37/5, lift 7.8)
## day <= 7
## month = aug
## poutcome = success
## -> class yes [0.846]
##
## Rule 24: (10/1, lift 7.7)
## marital = single
## housing = no
## day <= 9
## month = apr
## duration > 3
## -> class yes [0.833]
##
## Rule 25: (93/17, lift 7.5)
## duration > 3
## pdays > 385
## -> class yes [0.811]
##
## Rule 26: (7/1, lift 7.2)
## duration > 2
## campaign > 3
## previous > 5
## poutcome = success
## -> class yes [0.778]
##
## Rule 27: (47/10, lift 7.1)
## housing = no
## day > 9
## day <= 20
## month = feb
## duration > 1
## poutcome in {other, unknown}
## -> class yes [0.776]
##
## Rule 28: (20/4, lift 7.1)
## month = jun
## duration > 2
## pdays <= 385
## previous > 4
## poutcome in {failure, other}
## -> class yes [0.773]
##
## Rule 29: (434/99, lift 7.1)
## housing = no
## duration > 2
## poutcome = success
## -> class yes [0.771]
##
## Rule 30: (429/98, lift 7.1)
## duration > 3
## campaign <= 3
## poutcome = success
## -> class yes [0.770]
##
## Rule 31: (374/86, lift 7.1)
## month in {apr, aug, dec, jul, jun, mar, oct, sep}
## duration > 2
## campaign <= 3
## poutcome = success
## -> class yes [0.769]
##
## Rule 32: (154/36, lift 7.0)
## contact = cellular
## month in {mar, oct, sep}
## poutcome = success
## -> class yes [0.763]
##
## Rule 33: (187/50, lift 6.7)
## day > 20
## month in {apr, feb}
## duration > 3
## -> class yes [0.730]
##
## Rule 34: (5/1, lift 6.6)
## education = primary
## month in {mar, oct, sep}
## duration > 1
## previous > 4
## -> class yes [0.714]
##
## Rule 35: (195/63, lift 6.2)
## contact in {cellular, telephone}
## month = jun
## duration > 3
## pdays <= 385
## -> class yes [0.675]
##
## Rule 36: (608/235, lift 5.6)
## education in {secondary, tertiary}
## month in {dec, mar, oct, sep}
## duration > 2
## -> class yes [0.613]
##
## Rule 37: (289/120, lift 5.4)
## contact in {cellular, telephone}
## month = jun
## duration > 2
## pdays <= 385
## -> class yes [0.584]
##
## Rule 38: (204/93, lift 5.0)
## housing = no
## day > 9
## day <= 20
## month in {apr, feb}
## duration > 1
## -> class yes [0.544]
##
## Default class: no
##
##
## Evaluation on training data (26777 cases):
##
## Rules
## ----------------
## No Errors
##
## 38 2360( 8.8%) <<
##
##
## (a) (b) <-classified as
## ---- ----
## 23491 374 (a): class no
## 1986 926 (b): class yes
##
##
## Attribute usage:
##
## 97.21% poutcome
## 95.47% month
## 90.76% pdays
## 63.94% duration
## 63.44% housing
## 54.58% balance
## 51.17% previous
## 40.88% contact
## 31.44% day
## 17.03% education
## 12.87% campaign
## 4.09% marital
## 0.09% age
##
##
## Time: 0.8 secsL’arbre de decissió ha fet servir un total de 13 atributs per a generar les regles i ha generat un total de 38 regles.
set.seed(123)
model1 <- C50::C5.0(train.x, train.y, tree=TRUE)
summary(model1)##
## Call:
## C5.0.default(x = train.x, y = train.y, tree = TRUE)
##
##
## C5.0 [Release 2.07 GPL Edition] Tue Jun 16 22:30:13 2020
## -------------------------------
##
## Class specified by attribute `outcome'
##
## Read 26777 cases (17 attributes) from undefined.data
##
## Decision tree:
##
## poutcome = success:
## :...duration <= 2:
## : :...contact in {telephone,unknown}: no (13/1)
## : : contact = cellular:
## : : :...housing = yes: no (55/8)
## : : housing = no:
## : : :...month in {apr,dec,jan,jul,jun,may}: no (40/7)
## : : month in {mar,oct,sep}: yes (17/6)
## : : month = aug:
## : : :...day <= 7: yes (7/1)
## : : : day > 7: no (19/2)
## : : month = nov:
## : : :...balance <= 719: yes (4)
## : : : balance > 719: no (7/1)
## : : month = feb:
## : : :...pdays <= 92: yes (6)
## : : pdays > 92:
## : : :...age <= 3: no (5)
## : : age > 3: yes (4/1)
## : duration > 2:
## : :...housing = no: yes (434/99)
## : housing = yes:
## : :...campaign > 3:
## : :...previous <= 5: no (11/1)
## : : previous > 5: yes (5/1)
## : campaign <= 3:
## : :...month in {apr,aug,dec,jul,jun,mar,oct,sep}: yes (97/22)
## : month = jan:
## : :...pdays <= 155: yes (2)
## : : pdays > 155: no (3)
## : month = may:
## : :...duration <= 3: no (8/2)
## : : duration > 3: yes (52/22)
## : month = nov:
## : :...day <= 17: yes (10)
## : : day > 17:
## : : :...marital = divorced: yes (2)
## : : marital in {married,single}: no (10)
## : month = feb:
## : :...duration <= 3: no (4)
## : duration > 3:
## : :...previous > 3: yes (7)
## : previous <= 3:
## : :...day <= 4: no (4)
## : day > 4: yes (2)
## poutcome in {failure,other,unknown}:
## :...pdays > 385:
## :...duration <= 3: no (58/12)
## : duration > 3:
## : :...previous <= 4: yes (63/10)
## : previous > 4: no (10/3)
## pdays <= 385:
## :...month in {aug,jan,jul,jun,may,nov}:
## :...duration <= 2: no (8999/54)
## : duration > 2:
## : :...contact = unknown:
## : :...month in {jan,jul,jun,may}: no (4676/286)
## : : month in {aug,nov}:
## : : :...day <= 15: yes (7)
## : : day > 15: no (14/1)
## : contact in {cellular,telephone}:
## : :...month in {aug,jan,jul,may,nov}: no (7988/1109)
## : month = jun:
## : :...previous > 4: yes (20/4)
## : previous <= 4:
## : :...duration <= 3: no (75/25)
## : duration > 3: yes (143/52)
## month in {dec,mar,oct,sep}:
## :...duration <= 1: no (123/5)
## : duration > 1:
## : :...education = primary:
## : :...previous <= 4: no (53/16)
## : : previous > 4: yes (4/1)
## : education in {secondary,tertiary}:
## : :...duration > 2:
## : :...poutcome = failure: no (84/35)
## : : poutcome in {other,unknown}: yes (335/139)
## : duration <= 2:
## : :...previous > 2: no (18/1)
## : previous <= 2:
## : :...marital = divorced: yes (10/2)
## : marital in {married,single}: no (120/40)
## month in {apr,feb}:
## :...day > 20:
## :...duration <= 3: no (171/29)
## : duration > 3: yes (152/47)
## day <= 20:
## :...housing = yes: no (1849/125)
## housing = no:
## :...day > 9:
## :...poutcome = failure: no (25/1)
## : poutcome in {other,unknown}:
## : :...duration <= 1: no (17/2)
## : duration > 1:
## : :...month = apr: no (101/43)
## : month = feb: yes (47/10)
## day <= 9:
## :...month = feb: no (707/64)
## month = apr:
## :...loan = yes: yes (6/1)
## loan = no:
## :...marital in {divorced,married}: no (44/8)
## marital = single:
## :...duration > 3: yes (9/1)
## duration <= 3:
## :...balance <= 2174: no (17/2)
## balance > 2174: yes (4/1)
##
##
## Evaluation on training data (26777 cases):
##
## Decision Tree
## ----------------
## Size Errors
##
## 58 2303( 8.6%) <<
##
##
## (a) (b) <-classified as
## ---- ----
## 23445 420 (a): class no
## 1883 1029 (b): class yes
##
##
## Attribute usage:
##
## 100.00% poutcome
## 97.58% month
## 96.98% pdays
## 90.17% duration
## 48.92% contact
## 13.60% housing
## 12.04% day
## 2.33% education
## 2.04% previous
## 0.81% campaign
## 0.81% marital
## 0.30% loan
## 0.12% balance
## 0.03% age
##
##
## Time: 0.5 secsL’arbre té un tamany de 58 nodes amb una taxa d’error del 8.6%
predicted.model1 <- predict(model1, test.x, type="class")
print(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model1 == test.y) / length(predicted.model1)))## [1] "La precissió de l'arbre és del: 91.24 %"Tenim una precisió de model força bona, de més del 91% dels casos. Si per exemple, tinguéssim uns comercials que fessin campanyes telefòniques per a captar nous clients, podríem fer servir aquest model per a predir amb més d’un 90% de fiabilitat quin target atacar o quines decisions de comportament prendre.
Fent servir la matriu de confusió veiem la qualitat del model també:
## creem la matriu per avaluar la prova
matrix.conf <- table(Class=predicted.model1,Predicted=test.y)
percent.correct <- 100 * sum(diag(matrix.conf)) / sum(matrix.conf)
print(sprintf("L'error de classificació és: %.2f %%",100 - percent.correct))## [1] "L'error de classificació és: 8.76 %"mosaicplot(matrix.conf)
confusionMatrix(predicted.model1,test.y, dnn = c("Prediction"))## Confusion Matrix and Statistics
##
## NA
## Prediction no yes
## no 10066 812
## yes 193 405
##
## Accuracy : 0.9124
## 95% CI : (0.9071, 0.9175)
## No Information Rate : 0.894
## P-Value [Acc > NIR] : 2.267e-11
##
## Kappa : 0.4047
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.9812
## Specificity : 0.3328
## Pos Pred Value : 0.9254
## Neg Pred Value : 0.6773
## Prevalence : 0.8940
## Detection Rate : 0.8771
## Detection Prevalence : 0.9479
## Balanced Accuracy : 0.6570
##
## 'Positive' Class : no
## La veritat és que sense fer PCA tenim un valor de predicció molt bó amb un Kappa (aletorietat de la predicció) realment alt, cosa que es bó.
Fleiss et al. (2003) va afirmar que, a la majoria dels propòsits,
Altre cosa molt interessant:
Predim molt millor el ‘no’ que el ‘yes’. Aquesta part no m’agrada tant.
Preparem altre cop les dades per a avaluar l’arbre de decisió. Dividim en conjunt de prova i conjunt d’entrenament. Un el fem servir per a construir el model de l’arbre de decisió i el de prova per avaluar quina precisió obtenim d’aquest arbre. Fem servir la proporció 0.7 entrenament, 0.3 d’avaluació. La variable a avaluar en aquest cas és la columna ‘y’.
bank.data <- read.csv('bank-clean-discret.csv', sep = ',')
summary(bank.data)## age job marital education
## Min. :1.000 blue-collar:8507 divorced: 4494 primary : 5926
## 1st Qu.:2.000 management :7896 married :22762 secondary:20981
## Median :3.000 technician :6615 single :10997 tertiary :11346
## Mean :3.046 admin. :4565
## 3rd Qu.:4.000 services :3708
## Max. :5.000 retired :1492
## (Other) :5470
## default balance housing loan contact
## no :37485 Min. :-1884.0 no :16220 no :31548 cellular :24919
## yes: 768 1st Qu.: 43.0 yes:22033 yes: 6705 telephone: 2169
## Median : 338.0 unknown :11165
## Mean : 623.6
## 3rd Qu.: 957.0
## Max. : 3388.0
##
## day month duration campaign
## Min. : 1.00 may :12161 Min. :1.000 Min. : 1.000
## 1st Qu.: 8.00 jul : 6141 1st Qu.:2.000 1st Qu.: 1.000
## Median :16.00 aug : 5320 Median :3.000 Median : 2.000
## Mean :15.78 jun : 4333 Mean :3.003 Mean : 2.776
## 3rd Qu.:21.00 nov : 2940 3rd Qu.:4.000 3rd Qu.: 3.000
## Max. :31.00 apr : 2450 Max. :5.000 Max. :58.000
## (Other): 4908
## pdays previous poutcome y
## Min. : -1.0 Min. : 0.0000 failure: 4099 no :34124
## 1st Qu.: -1.0 1st Qu.: 0.0000 other : 1551 yes: 4129
## Median : -1.0 Median : 0.0000 success: 1156
## Mean : 40.3 Mean : 0.5659 unknown:31447
## 3rd Qu.: -1.0 3rd Qu.: 0.0000
## Max. :871.0 Max. :275.0000
## # separem el 70% de conjunt d'entrenament del 30% per a avaluar
## obetnim primer una mostra del total de forma aleatoria (hem guardat físicament aquest random per poder reproduïr el mateix i el carreguem de disc)
random.data <- read.csv('random-data.csv', sep = ',')
random.data <- unlist(random.data, use.names=FALSE)
data.train <- bank.data[random.data,]
data.test <- bank.data[-random.data,]
# en les dades d'entrenament, separem les dades per a generar l'arbre de la dada objectiu (target)
train.x <- data.train[,1:16]
train.y <- data.train[,17]
# en les dades de test, separem les dades per a generar l'arbre de la dada objectiu (target)
test.x <- data.test[,1:16]
test.y <- data.test[,17]L’executem de forma que ens generi regles per a fer prediccions. Ara executem l’algoritme amb l’opció Winnowed. Aquesta opció no farà servir aquelles dimensions que no trobi que siguin rellevants, es a dir, està fent PCA de forma transparent i automàtica.
# amb l'atribut winnow eliminem dimensions que no tenen impacte en les classes de sortida
model2 <- C50::C5.0(train.x, train.y, rules=TRUE, size=4, control = C5.0Control(winnow = TRUE))
summary(model2)##
## Call:
## C5.0.default(x = train.x, y = train.y, rules = TRUE, control
## = C5.0Control(winnow = TRUE), size = 4)
##
##
## C5.0 [Release 2.07 GPL Edition] Tue Jun 16 22:30:16 2020
## -------------------------------
##
## Class specified by attribute `outcome'
##
## Read 26777 cases (17 attributes) from undefined.data
##
## 7 attributes winnowed
## Estimated importance of remaining attributes:
##
## 21% contact
## 6% poutcome
## 3% month
## 2% duration
## 1% housing
## 1% day
## 1% loan
## 1% pdays
## 1% balance
##
## Rules:
##
## Rule 1: (152, lift 1.1)
## housing = yes
## month = jan
## pdays > 155
## -> class no [0.994]
##
## Rule 2: (4295/27, lift 1.1)
## housing = yes
## day <= 19
## month = may
## duration <= 4
## -> class no [0.993]
##
## Rule 3: (10681/228, lift 1.1)
## duration <= 2
## -> class no [0.979]
##
## Rule 4: (7737/302, lift 1.1)
## contact = unknown
## month in {aug, jan, jul, jun, may, nov}
## poutcome in {failure, other, unknown}
## -> class no [0.961]
##
## Rule 5: (930/45, lift 1.1)
## housing = yes
## day > 17
## month = nov
## -> class no [0.951]
##
## Rule 6: (256/13, lift 1.1)
## housing = yes
## loan = yes
## month = jun
## -> class no [0.946]
##
## Rule 7: (606/35, lift 1.1)
## housing = yes
## day <= 7
## month = feb
## -> class no [0.941]
##
## Rule 8: (1183/79, lift 1.0)
## housing = yes
## day <= 20
## month = apr
## -> class no [0.932]
##
## Rule 9: (25818/2325, lift 1.0)
## pdays <= 385
## poutcome in {failure, other, unknown}
## -> class no [0.910]
##
## Rule 10: (9/1, lift 7.5)
## contact = cellular
## day <= 5
## month = dec
## duration > 2
## poutcome in {failure, other, unknown}
## -> class yes [0.818]
##
## Rule 11: (93/17, lift 7.5)
## duration > 3
## pdays > 385
## -> class yes [0.811]
##
## Rule 12: (47/10, lift 7.1)
## housing = no
## day > 9
## day <= 20
## month = feb
## duration > 1
## poutcome in {other, unknown}
## -> class yes [0.776]
##
## Rule 13: (13/3, lift 6.7)
## balance > 326
## balance <= 1059
## contact = cellular
## month = dec
## duration > 2
## poutcome in {other, unknown}
## -> class yes [0.733]
##
## Rule 14: (135/37, lift 6.6)
## month = sep
## duration > 3
## -> class yes [0.723]
##
## Rule 15: (88/24, lift 6.6)
## loan = no
## month = mar
## duration > 2
## poutcome in {other, unknown}
## -> class yes [0.722]
##
## Rule 16: (5/1, lift 6.6)
## day > 3
## day <= 4
## month in {jun, nov}
## duration > 2
## poutcome = failure
## -> class yes [0.714]
##
## Rule 17: (83/24, lift 6.5)
## contact in {cellular, telephone}
## day <= 3
## month in {jul, jun, nov}
## duration > 3
## poutcome = unknown
## -> class yes [0.706]
##
## Rule 18: (13/4, lift 6.1)
## balance > 1784
## month = dec
## duration > 1
## poutcome in {failure, other, unknown}
## -> class yes [0.667]
##
## Rule 19: (141/48, lift 6.0)
## day > 15
## month = oct
## duration > 2
## -> class yes [0.657]
##
## Rule 20: (828/309, lift 5.8)
## poutcome = success
## -> class yes [0.627]
##
## Rule 21: (657/260, lift 5.6)
## month in {dec, mar, oct, sep}
## duration > 2
## -> class yes [0.604]
##
## Rule 22: (161/68, lift 5.3)
## pdays > 385
## -> class yes [0.577]
##
## Rule 23: (204/93, lift 5.0)
## housing = no
## day > 9
## day <= 20
## month in {apr, feb}
## duration > 1
## -> class yes [0.544]
##
## Rule 24: (230/105, lift 5.0)
## contact in {cellular, telephone}
## day <= 4
## month in {jul, jun, nov}
## duration > 2
## -> class yes [0.543]
##
## Default class: no
##
##
## Evaluation on training data (26777 cases):
##
## Rules
## ----------------
## No Errors
##
## 24 2385( 8.9%) <<
##
##
## (a) (b) <-classified as
## ---- ----
## 23573 292 (a): class no
## 2093 819 (b): class yes
##
##
## Attribute usage:
##
## 99.55% poutcome
## 97.03% pdays
## 52.73% duration
## 52.62% month
## 29.83% contact
## 28.48% housing
## 28.38% day
## 1.28% loan
## 0.10% balance
##
##
## Time: 0.7 secspredicted.model2 <- predict(model2, test.x, type="class")
cat(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model2 == test.y) / length(predicted.model2)))## La precissió de l'arbre és del: 91.06 %Anem a veure els resultats:
S’han llegit 26777 casos amb 9 atributs. Amb el winnowed s’han descartat 7 atributs per a construir les regles.
Classe per defecte -> ‘no’ (no contracta crèdit bancari)
S’han generat 24 regles de classificació
Ha fet servir els atributs:
L’atribut amb més pes i que és l’arrel de l’arbre es poutcome.
L’arbre ha classificat malament 24 de 2385 cassos (8.9% d’errors):
Tots els lift, aleatorietat de la predicció estan per sobre d’1, per tant les regles són prou vàlides.
L’arbre ha generat 24 regles de decisió. Podem observar que la fiabilitat de les primeres regles son totes les que estan orientades a la resposta ‘no’. Ens es molt més fàcil predir el ‘no’ que el ‘yes’. La primera regla amb ‘yes’, és la 10 amb un confidence de 0.818.
set.seed(123)
model2 <- C50::C5.0(train.x, train.y, control = C5.0Control(winnow = TRUE))
plot(model2)
summary(model2)##
## Call:
## C5.0.default(x = train.x, y = train.y, control = C5.0Control(winnow = TRUE))
##
##
## C5.0 [Release 2.07 GPL Edition] Tue Jun 16 22:30:19 2020
## -------------------------------
##
## Class specified by attribute `outcome'
##
## Read 26777 cases (17 attributes) from undefined.data
##
## 7 attributes winnowed
## Estimated importance of remaining attributes:
##
## 21% contact
## 6% poutcome
## 3% month
## 2% duration
## 1% housing
## 1% day
## 1% loan
## 1% pdays
## 1% balance
##
## Decision tree:
##
## poutcome = success:
## :...duration <= 2: no (177/49)
## : duration > 2:
## : :...housing = no: yes (434/99)
## : housing = yes:
## : :...month in {aug,dec,jul,mar,oct,sep}: yes (67/11)
## : month = apr:
## : :...day <= 20: no (14/4)
## : : day > 20: yes (12/2)
## : month = jan:
## : :...pdays <= 155: yes (2)
## : : pdays > 155: no (3)
## : month = jun:
## : :...loan = no: yes (9/1)
## : : loan = yes: no (2)
## : month = nov:
## : :...day <= 17: yes (10)
## : : day > 17: no (12/2)
## : month = feb:
## : :...duration <= 3: no (4)
## : : duration > 3:
## : : :...day <= 7: no (8/3)
## : : day > 7: yes (6)
## : month = may:
## : :...duration > 4: yes (34/12)
## : duration <= 4:
## : :...day <= 19: no (30/8)
## : day > 19: yes (4)
## poutcome in {failure,other,unknown}:
## :...pdays > 385:
## :...duration <= 3: no (58/12)
## : duration > 3: yes (73/17)
## pdays <= 385:
## :...month in {apr,feb}:
## :...day > 20:
## : :...duration <= 3: no (171/29)
## : : duration > 3: yes (152/47)
## : day <= 20:
## : :...housing = yes: no (1849/125)
## : housing = no:
## : :...day <= 9:
## : :...month = feb: no (707/64)
## : : month = apr:
## : : :...loan = no: no (74/21)
## : : loan = yes: yes (6/1)
## : day > 9:
## : :...poutcome = failure: no (25/1)
## : poutcome in {other,unknown}:
## : :...duration <= 1: no (17/2)
## : duration > 1:
## : :...month = apr: no (101/43)
## : month = feb: yes (47/10)
## month in {aug,jan,jul,jun,may,nov}:
## :...day > 4: no (20606/1395)
## : day <= 4:
## : :...contact = unknown: no (775/33)
## : contact in {cellular,telephone}:
## : :...duration <= 2: no (183/12)
## : duration > 2:
## : :...month in {aug,jan,may}: no (155/44)
## : month in {jul,jun,nov}:
## : :...duration <= 3:
## : :...day <= 2: no (42/11)
## : : day > 2: yes (29/13)
## : duration > 3:
## : :...poutcome = other: no (7/2)
## : poutcome = failure:
## : :...day <= 3: no (24/10)
## : : day > 3: yes (4/1)
## : poutcome = unknown:
## : :...day <= 3: yes (83/24)
## : day > 3: no (14/4)
## month in {dec,mar,oct,sep}:
## :...duration <= 1: no (123/5)
## duration > 1:
## :...duration <= 2: no (163/50)
## duration > 2:
## :...poutcome = failure: no (93/39)
## poutcome in {other,unknown}:
## :...month = mar:
## :...loan = no: yes (87/24)
## : loan = yes: no (6/2)
## month = oct:
## :...day <= 15: no (44/11)
## : day > 15: yes (85/30)
## month = sep:
## :...duration <= 3: no (24/8)
## : duration > 3: yes (71/25)
## month = dec:
## :...contact in {telephone,unknown}: no (4)
## contact = cellular:
## :...day <= 5: yes (7/1)
## day > 5:
## :...balance > 1784: yes (8/3)
## balance <= 1784:
## :...balance > 1059: no (8/1)
## balance <= 1059:
## :...balance <= 326: no (13/5)
## balance > 326: yes (11/3)
##
##
## Evaluation on training data (26777 cases):
##
## Decision Tree
## ----------------
## Size Errors
##
## 55 2319( 8.7%) <<
##
##
## (a) (b) <-classified as
## ---- ----
## 23541 324 (a): class no
## 1995 917 (b): class yes
##
##
## Attribute usage:
##
## 100.00% poutcome
## 97.23% month
## 96.93% pdays
## 94.64% day
## 12.99% housing
## 10.21% duration
## 5.11% contact
## 0.69% loan
## 0.15% balance
##
##
## Time: 0.5 secsHe generat un arbre de 55 nodes. En l’arbre podem veure que el primer nivell es realitza amb l’atribut poutcome. És el que aporta particions més homogènies.
predicted.model2 <- predict(model2, test.x, type="class")
print(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model2 == test.y) / length(predicted.model2)))## [1] "La precissió de l'arbre és del: 91.30 %"Tenin una precissió de model força bona, de mes del 91% dels casos, però està una mica per sota que sense fer winnowed.
Fent servir la matriu de confusió veiem la qualitat del model també:
## creem la matriu per avaluar la prova
matrix.conf <- table(Class=predicted.model2,Predicted=test.y)
percent.correct <- 100 * sum(diag(matrix.conf)) / sum(matrix.conf)
print(sprintf("L'error de classificació és: %.2f %%",100 - percent.correct))## [1] "L'error de classificació és: 8.70 %"mosaicplot(matrix.conf)
Encara que la predicció sigui realment bona hi ha una cosa que no m’acaba d’agradar i és que en la predicció del ‘yes’ genera molt mes errors. Quan tenim un cas que ens diu que pot contractar un disposit bancari, té força errors. En canvi, el ‘no’, és molt precís. De totes maneres, és lleugerament millor que sense ‘winnow’.
confusionMatrix(predicted.model2,test.y, dnn = c("Prediction"))## Confusion Matrix and Statistics
##
## NA
## Prediction no yes
## no 10116 855
## yes 143 362
##
## Accuracy : 0.913
## 95% CI : (0.9077, 0.9181)
## No Information Rate : 0.894
## P-Value [Acc > NIR] : 4.875e-12
##
## Kappa : 0.382
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.9861
## Specificity : 0.2975
## Pos Pred Value : 0.9221
## Neg Pred Value : 0.7168
## Prevalence : 0.8940
## Detection Rate : 0.8815
## Detection Prevalence : 0.9560
## Balanced Accuracy : 0.6418
##
## 'Positive' Class : no
## Amb les estadístiques veiem que el ‘no’ té un 0.9221 i en canvi el ‘yes’ 0.7168 El valor de Kappa que tenim no és molt bo. Aquest valor és una altra forma d’interpretar l’aleatorietat explicat en el punt anterior. Ens ha baixat una mica fent servir menys atributs.
Havíem reservat un model per a provar si millorava o empitjorava el fet d’eliminar outliers i alguns camps nulls.
bank.data <- read.csv('bank-discret.csv', sep = ',')
# separem el 70% de conjunt d'entrenament del 30% per a avaluar
random.data <- sample(1:nrow(bank.data), 0.70 * nrow(bank.data))
data.train <- bank.data[random.data,]
data.test <- bank.data[-random.data,]
# en les dades d'entrenament, separem les dades per a generar l'arbre de la dada objectiu (target)
train.x <- data.train[,1:16]
train.y <- data.train[,17]
# en les dades de test, separem les dades per a generar l'arbre de la dada objectiu (target)
test.x <- data.test[,1:16]
test.y <- data.test[,17]
model3 <- C50::C5.0(train.x, train.y, rules=TRUE, size=10,control = C5.0Control(winnow = TRUE))
predicted.model3 <- predict(model3, test.x, type="class")
matrix.conf <- table(Class=predicted.model3,Predicted=test.y)
mosaicplot(matrix.conf)
confusionMatrix(predicted.model3,test.y, dnn = c("Prediction"))## Confusion Matrix and Statistics
##
## NA
## Prediction no yes
## no 11716 1087
## yes 263 498
##
## Accuracy : 0.9005
## 95% CI : (0.8953, 0.9055)
## No Information Rate : 0.8831
## P-Value [Acc > NIR] : 7.382e-11
##
## Kappa : 0.3773
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.9780
## Specificity : 0.3142
## Pos Pred Value : 0.9151
## Neg Pred Value : 0.6544
## Prevalence : 0.8831
## Detection Rate : 0.8638
## Detection Prevalence : 0.9439
## Balanced Accuracy : 0.6461
##
## 'Positive' Class : no
## La predicció en quant a accuracy és pitjor, el valor kappa (aletorietat) també és pitjor. La predicció del ‘yes’ ens cau a ‘0.6544’. Están ben netejades les dades perquè ens perjudica al classificar tenir valors nulls i outliers.
Anem a reduïr encara més la dimensionalitat del dataset. Hi ha alguns atributs que potser ens fan embrutar l’arbre de decissió amb regles no necessaries
Primer veiem la matriu de correlació per veure quins atributs tenen impacte amb la variable de sortida target.
bank.data <- read.csv('bank-clean-discret.csv', sep = ',')
bank.data.tmp <- bank.data
for (colname in colnames(bank.data.tmp)) {
bank.data.tmp[colname] <- lapply(bank.data.tmp[colname],as.integer)
}
corr.mat <- cor(bank.data.tmp)
# visualize it
library(corrplot)## Warning: package 'corrplot' was built under R version 3.6.3## corrplot 0.84 loadedcol <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))
corrplot(corr.mat, method="color", col=col(200),
addCoef.col = "black",
tl.col="black", tl.srt=45,
insig = "blank",
number.cex=0.5)
Veiem que els atributs que tenen més impacte en quant a la matriu de correlació serien (agafem > 0.07): education, balance, housing, contact, duration, pdays, previous
En l’estudi de PCA el que es vol es trobar un nombre de dimensions menor de la mostra que expliquin el mateix que la mostra total. A nivell matemàtic funciona amb el concepte de ‘eigenvectors’ i ‘eigenvalues’, que són valors que simplifiquen una matriu dient el mateix que la matriu original: Eigenvector and Eigenvalue
El que fem és trobar la variànça total del conjunt de dades i intentarem reduir la dimensionalitat mantenint aquesta varinça al màxim. Per exemple, si tinguesim amb 10 dimensions un > 95% de la variança ja ens donariem per satisfets. L’elecció es realitza de manera que la primera component principal sigui la que major variància reculli; la segona ha de recollir la màxima variabilitat no recollida per la primera, i així successivament, triant un nombre que reculli un percentatge suficient de variància total.
Dividirem l’estudi en les dades quantitatives (PCA, Principal Component Analysis) i mixtes (FAMD, Factor Analysis of Mixed Data) per a dades qualitatives i quantitatives. Info
bank.pca <- read.csv('bank-clean-discret.csv', sep = ',')
quantitative <- bank.pca[,c(1,6,10,12:15)]
# executem l'estudi PCA
pca <- prcomp(quantitative, scale. = TRUE)
# la funció mostra la variança de cada component
fviz_eig(pca, choice = c("variance"), addlabels = TRUE)
# relació entre les dimensions i atributs.
fviz_contrib(pca, choice = "var", axes = 1:4)
# pve acumulada
PVE <- 100*pca$sdev^2/sum(pca$sdev^2)
plot(cumsum(PVE), type = "o",
ylab = "PVE acumulada",
xlab = "Componente principal",
col = "blue")
Totes les variables continues tenen una bona aportació encara que podriem reduïr el nombre ja que amb 4 dimensions estem aportant molta informació de la mostra.
En la tercera gràfica es pot veure que amb 5 dimensions tenim una variança acumulada de quasi un 80%. Podríem eliminar day, duration i campaign
## obetnim primer una mostra del total de forma aleatoria (hem guardat físicament aquest random per poder reproduïr el mateix i el carreguem de disc)
bank.data <- read.csv('bank-clean-discret.csv', sep = ',')
summary(bank.data)## age job marital education
## Min. :1.000 blue-collar:8507 divorced: 4494 primary : 5926
## 1st Qu.:2.000 management :7896 married :22762 secondary:20981
## Median :3.000 technician :6615 single :10997 tertiary :11346
## Mean :3.046 admin. :4565
## 3rd Qu.:4.000 services :3708
## Max. :5.000 retired :1492
## (Other) :5470
## default balance housing loan contact
## no :37485 Min. :-1884.0 no :16220 no :31548 cellular :24919
## yes: 768 1st Qu.: 43.0 yes:22033 yes: 6705 telephone: 2169
## Median : 338.0 unknown :11165
## Mean : 623.6
## 3rd Qu.: 957.0
## Max. : 3388.0
##
## day month duration campaign
## Min. : 1.00 may :12161 Min. :1.000 Min. : 1.000
## 1st Qu.: 8.00 jul : 6141 1st Qu.:2.000 1st Qu.: 1.000
## Median :16.00 aug : 5320 Median :3.000 Median : 2.000
## Mean :15.78 jun : 4333 Mean :3.003 Mean : 2.776
## 3rd Qu.:21.00 nov : 2940 3rd Qu.:4.000 3rd Qu.: 3.000
## Max. :31.00 apr : 2450 Max. :5.000 Max. :58.000
## (Other): 4908
## pdays previous poutcome y
## Min. : -1.0 Min. : 0.0000 failure: 4099 no :34124
## 1st Qu.: -1.0 1st Qu.: 0.0000 other : 1551 yes: 4129
## Median : -1.0 Median : 0.0000 success: 1156
## Mean : 40.3 Mean : 0.5659 unknown:31447
## 3rd Qu.: -1.0 3rd Qu.: 0.0000
## Max. :871.0 Max. :275.0000
## bank.data.subset <- bank.data[,c("balance","age","previous","pdays","day","y")]
summary(bank.data.subset)## balance age previous pdays
## Min. :-1884.0 Min. :1.000 Min. : 0.0000 Min. : -1.0
## 1st Qu.: 43.0 1st Qu.:2.000 1st Qu.: 0.0000 1st Qu.: -1.0
## Median : 338.0 Median :3.000 Median : 0.0000 Median : -1.0
## Mean : 623.6 Mean :3.046 Mean : 0.5659 Mean : 40.3
## 3rd Qu.: 957.0 3rd Qu.:4.000 3rd Qu.: 0.0000 3rd Qu.: -1.0
## Max. : 3388.0 Max. :5.000 Max. :275.0000 Max. :871.0
## day y
## Min. : 1.00 no :34124
## 1st Qu.: 8.00 yes: 4129
## Median :16.00
## Mean :15.78
## 3rd Qu.:21.00
## Max. :31.00random.data <- read.csv('random-data.csv', sep = ',')
random.data <- unlist(random.data, use.names=FALSE)
data.train <- bank.data.subset[random.data,]
data.test <- bank.data.subset[-random.data,]
# en les dades d'entrenament, separem les dades per a generar l'arbre de la dada objectiu (target)
train.x <- data.train[,1:5]
train.y <- data.train[,6]
# en les dades de test, separem les dades per a generar l'arbre de la dada objectiu (target)
test.x <- data.test[,1:5]
test.y <- data.test[,6]
model31 <- C50::C5.0(train.x, train.y)
predicted.model31 <- predict(model31, test.x, type="class")
print(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model31 == test.y) / length(predicted.model31)))## [1] "La precissió de l'arbre és del: 89.85 %"matrix.conf <- table(Class=predicted.model31,Predicted=test.y)
percent.correct <- 100 * sum(diag(matrix.conf)) / sum(matrix.conf)
print(sprintf("L'error de classificació és: %.2f %%",100 - percent.correct))## [1] "L'error de classificació és: 10.15 %"mosaicplot(matrix.conf)
confusionMatrix(predicted.model31,test.y, dnn = c("Prediction"))## Confusion Matrix and Statistics
##
## NA
## Prediction no yes
## no 10150 1056
## yes 109 161
##
## Accuracy : 0.8985
## 95% CI : (0.8928, 0.904)
## No Information Rate : 0.894
## P-Value [Acc > NIR] : 0.05853
##
## Kappa : 0.1852
##
## Mcnemar's Test P-Value : < 2e-16
##
## Sensitivity : 0.9894
## Specificity : 0.1323
## Pos Pred Value : 0.9058
## Neg Pred Value : 0.5963
## Prevalence : 0.8940
## Detection Rate : 0.8845
## Detection Prevalence : 0.9765
## Balanced Accuracy : 0.5608
##
## 'Positive' Class : no
## La predicció no és gens bona i tenim un Kappa molt baix. Realment es troben a faltar les variables categòr
Com que el nostre dataset té dades categòriques i quantitatives, anem a fer l’estudi de dades mixtes amb la funció FAMD.
mixed <- read.csv('bank-clean-discret.csv', sep = ',')
res.famd <- FAMD(mixed, sup.var = 17, graph = FALSE, ncp = 40)
# relació entre les dimensions i atributs.
fviz_contrib(res.famd, choice = "var", axes = 1:40)
plot(res.famd$eig[,3], type = "o",
ylab = "PVE acumulada",
xlab = "Componente principal",
col = "blue")
Veiem en l’estudi que tenim 5 o 6 components que aporten moltíssim i la resta tenen un impacte similar. Anem a agafar un subconjunt del dataset que tenim guardat només amb els atributs: month, job, education, poutcome i contact.
## obetnim primer una mostra del total de forma aleatoria (hem guardat físicament aquest random per poder reproduïr el mateix i el carreguem de disc)
bank.data <- read.csv('bank-clean-discret.csv', sep = ',')
summary(bank.data)## age job marital education
## Min. :1.000 blue-collar:8507 divorced: 4494 primary : 5926
## 1st Qu.:2.000 management :7896 married :22762 secondary:20981
## Median :3.000 technician :6615 single :10997 tertiary :11346
## Mean :3.046 admin. :4565
## 3rd Qu.:4.000 services :3708
## Max. :5.000 retired :1492
## (Other) :5470
## default balance housing loan contact
## no :37485 Min. :-1884.0 no :16220 no :31548 cellular :24919
## yes: 768 1st Qu.: 43.0 yes:22033 yes: 6705 telephone: 2169
## Median : 338.0 unknown :11165
## Mean : 623.6
## 3rd Qu.: 957.0
## Max. : 3388.0
##
## day month duration campaign
## Min. : 1.00 may :12161 Min. :1.000 Min. : 1.000
## 1st Qu.: 8.00 jul : 6141 1st Qu.:2.000 1st Qu.: 1.000
## Median :16.00 aug : 5320 Median :3.000 Median : 2.000
## Mean :15.78 jun : 4333 Mean :3.003 Mean : 2.776
## 3rd Qu.:21.00 nov : 2940 3rd Qu.:4.000 3rd Qu.: 3.000
## Max. :31.00 apr : 2450 Max. :5.000 Max. :58.000
## (Other): 4908
## pdays previous poutcome y
## Min. : -1.0 Min. : 0.0000 failure: 4099 no :34124
## 1st Qu.: -1.0 1st Qu.: 0.0000 other : 1551 yes: 4129
## Median : -1.0 Median : 0.0000 success: 1156
## Mean : 40.3 Mean : 0.5659 unknown:31447
## 3rd Qu.: -1.0 3rd Qu.: 0.0000
## Max. :871.0 Max. :275.0000
## bank.data.subset <- bank.data[,c("month","job","education","poutcome","contact","y")]
summary(bank.data.subset)## month job education poutcome
## may :12161 blue-collar:8507 primary : 5926 failure: 4099
## jul : 6141 management :7896 secondary:20981 other : 1551
## aug : 5320 technician :6615 tertiary :11346 success: 1156
## jun : 4333 admin. :4565 unknown:31447
## nov : 2940 services :3708
## apr : 2450 retired :1492
## (Other): 4908 (Other) :5470
## contact y
## cellular :24919 no :34124
## telephone: 2169 yes: 4129
## unknown :11165
##
##
##
## random.data <- read.csv('random-data.csv', sep = ',')
random.data <- unlist(random.data, use.names=FALSE)
data.train <- bank.data.subset[random.data,]
data.test <- bank.data.subset[-random.data,]
# en les dades d'entrenament, separem les dades per a generar l'arbre de la dada objectiu (target)
train.x <- data.train[,1:5]
train.y <- data.train[,6]
# en les dades de test, separem les dades per a generar l'arbre de la dada objectiu (target)
test.x <- data.test[,1:5]
test.y <- data.test[,6]
model4 <- C50::C5.0(train.x, train.y, rules=TRUE)
summary(model4)##
## Call:
## C5.0.default(x = train.x, y = train.y, rules = TRUE)
##
##
## C5.0 [Release 2.07 GPL Edition] Tue Jun 16 22:30:50 2020
## -------------------------------
##
## Class specified by attribute `outcome'
##
## Read 26777 cases (6 attributes) from undefined.data
##
## Rules:
##
## Rule 1: (25949/2393, lift 1.0)
## poutcome in {failure, other, unknown}
## -> class no [0.908]
##
## Rule 2: (828/309, lift 5.8)
## poutcome = success
## -> class yes [0.627]
##
## Default class: no
##
##
## Evaluation on training data (26777 cases):
##
## Rules
## ----------------
## No Errors
##
## 2 2702(10.1%) <<
##
##
## (a) (b) <-classified as
## ---- ----
## 23556 309 (a): class no
## 2393 519 (b): class yes
##
##
## Attribute usage:
##
## 100.00% poutcome
##
##
## Time: 0.1 secspredicted.model4 <- predict(model4, test.x, type="class")
print(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model4 == test.y) / length(predicted.model4)))## [1] "La precissió de l'arbre és del: 90.39 %"Obtenim un arbre molt més simple amb unes regles també més simples però el valor de la predicció ens ha baixar una mica.
predicted.model4 <- predict(model4, test.x, type="class")
print(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model4 == test.y) / length(predicted.model4)))## [1] "La precissió de l'arbre és del: 90.39 %"Amb una predicció lleugerament inferior al model inicial, però no moltíssim.
matrix.conf <- table(Class=predicted.model4,Predicted=test.y)
percent.correct <- 100 * sum(diag(matrix.conf)) / sum(matrix.conf)
print(sprintf("L'error de classificació és: %.2f %%",100 - percent.correct))## [1] "L'error de classificació és: 9.61 %"mosaicplot(matrix.conf)
confusionMatrix(predicted.model4,test.y, dnn = c("Prediction"))## Confusion Matrix and Statistics
##
## NA
## Prediction no yes
## no 10152 996
## yes 107 221
##
## Accuracy : 0.9039
## 95% CI : (0.8983, 0.9092)
## No Information Rate : 0.894
## P-Value [Acc > NIR] : 0.000245
##
## Kappa : 0.2524
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.9896
## Specificity : 0.1816
## Pos Pred Value : 0.9107
## Neg Pred Value : 0.6738
## Prevalence : 0.8940
## Detection Rate : 0.8846
## Detection Prevalence : 0.9714
## Balanced Accuracy : 0.5856
##
## 'Positive' Class : no
## El que si que ens ha perjudicat és en el valor de Kappa, on sembla que l’aletorietat és bastant pitjor
model4 <- C50::C5.0(train.x, train.y, rules=FALSE)
plot(model4)
Hem fet 4 estudis:
Sense PCA:
Amb PCA mitjançant winnowed
Amb PCA manual mitjançant estudi dels PC:
Sense eliminar outliers ni valors nulls:
Conclusions: - Podem dir que la classificació amb PCA(winnowed) i sense PCA ens ha donat valors realment molt similars. - Amb PCA (winnowed), hem millorat en la predicció general a causa que tenim millor predicció del ‘yes’. - Sense PCA hem obtingut bons valors i el millor valor de Kappa. - La millora en la predicció aplicant PCA(winnowed) és molt petita. Això pot venir donat per el fet que necessitem moltes dimensions per a tenir una variança acumulada per sobre de 95%. No podem desfer-nos de moltes varibles.
Un Random Forest és un conjunt d’arbres de decisió combinats amb bagging. A l’usar bagging, el que en realitat està passant, és que diferents arbres veuen diferents porcions de les dades. Cap arbre veu totes les dades d’entrenament. Això fa que cada arbre s’entreni amb diferents mostres de dades per a un mateix problema. D’aquesta manera, a l’combinar els seus resultats, uns errors es compensen amb altres i tenim una predicció que generalitza millor.
bank.data <- read.csv('bank-clean-discret.csv', sep = ',')
random.data <- read.csv('random-data.csv', sep = ',')
random.data <- unlist(random.data, use.names=FALSE)
train <- bank.data[random.data,]
test <- bank.data[-random.data,]
test.x <- test[,1:16]
test.y <- test[,17]
train_rf <- randomForest(y~.,data = train, ntree = 500)
predicted.model5 <- predict(train_rf,test.x)
print(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model5 == test.y) / length(predicted.model5)))## [1] "La precissió de l'arbre és del: 91.08 %"matrix.conf <- table(Class=predicted.model5,Predicted=test.y)
percent.correct <- 100 * sum(diag(matrix.conf)) / sum(matrix.conf)
print(sprintf("L'error de classificació és: %.2f %%",100 - percent.correct))## [1] "L'error de classificació és: 8.92 %"mosaicplot(matrix.conf)
confusionMatrix(predicted.model5,test.y, dnn = c("Prediction"))## Confusion Matrix and Statistics
##
## NA
## Prediction no yes
## no 10033 798
## yes 226 419
##
## Accuracy : 0.9108
## 95% CI : (0.9054, 0.9159)
## No Information Rate : 0.894
## P-Value [Acc > NIR] : 1.13e-09
##
## Kappa : 0.4064
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.9780
## Specificity : 0.3443
## Pos Pred Value : 0.9263
## Neg Pred Value : 0.6496
## Prevalence : 0.8940
## Detection Rate : 0.8743
## Detection Prevalence : 0.9438
## Balanced Accuracy : 0.6611
##
## 'Positive' Class : no
## Random Forest:
Tenim bona predicció i Kappa però la predicció del ‘no’ és força baixa (0.6496)
A diferència del ID3, aquest tipus d’arbre:
Mesura d’homogeneïtat i criteri d’aturada en el procés de partició i divisió de l’arbre a partir de l’índex de Gini, encara que hi ha variants que n’escullen d’altres. Aquesta mesura estableix el millor separador com el que redueix la diversitat de les diferents particions obtingudes (subarbres).
La poda es realitza de la següent forma:
Són binaris; a cada node hi ha un punt de tall (per un procediment semblant al que s’ha explicat per a trobar punts de tall en la discretització d’atributs continus) que separa en dos el conjunt d’observacions.
L’algorisme CART pot treballar amb atributs continus (tot i que les modificacions de ID3 també ho poden fer com hem pogut veure en l’exemple anterior C4.5).
Pot fer tant classificació com regressió: en el primer cas, la variable per a predir ha de ser categòrica amb un valor per a cada classe possible.
bank.data <- read.csv('bank-clean-discret.csv', sep = ',')
random.data <- read.csv('random-data.csv', sep = ',')
random.data <- unlist(random.data, use.names=FALSE)
train <- bank.data[random.data,]
test <- bank.data[-random.data,]
test.x <- test[,1:16]
test.y <- test[,17]
heart.tree <- rpart(y ~ .,method="class", data=train)
summary(heart.tree)## Call:
## rpart(formula = y ~ ., data = train, method = "class")
## n= 26777
##
## CP nsplit rel error xerror xstd
## 1 0.02403846 0 1.0000000 1.0000000 0.01749460
## 2 0.01236264 4 0.9007555 0.9007555 0.01670403
## 3 0.01000000 6 0.8760302 0.8839286 0.01656405
##
## Variable importance
## duration poutcome month contact
## 53 36 6 5
##
## Node number 1: 26777 observations, complexity param=0.02403846
## predicted class=no expected loss=0.10875 P(node) =1
## class counts: 23865 2912
## probabilities: 0.891 0.109
## left son=2 (21471 obs) right son=3 (5306 obs)
## Primary splits:
## duration < 4.5 to the left, improve=597.0337, (0 missing)
## poutcome splits as LLRL, improve=458.6320, (0 missing)
## month splits as LLRLLLLRLLRR, improve=270.8700, (0 missing)
## pdays < 8.5 to the left, improve=131.0225, (0 missing)
## previous < 0.5 to the left, improve=129.9314, (0 missing)
## Surrogate splits:
## balance < 3386.5 to the left, agree=0.802, adj=0, (0 split)
## previous < 53 to the left, agree=0.802, adj=0, (0 split)
##
## Node number 2: 21471 observations, complexity param=0.02403846
## predicted class=no expected loss=0.05626193 P(node) =0.8018449
## class counts: 20263 1208
## probabilities: 0.944 0.056
## left son=4 (20854 obs) right son=5 (617 obs)
## Primary splits:
## poutcome splits as LLRL, improve=342.36130, (0 missing)
## month splits as LLRLLLLRLLRR, improve=216.18880, (0 missing)
## pdays < 0 to the left, improve= 98.61017, (0 missing)
## previous < 0.5 to the left, improve= 98.61017, (0 missing)
## housing splits as RL, improve= 53.97091, (0 missing)
##
## Node number 3: 5306 observations, complexity param=0.02403846
## predicted class=no expected loss=0.3211459 P(node) =0.1981551
## class counts: 3602 1704
## probabilities: 0.679 0.321
## left son=6 (5095 obs) right son=7 (211 obs)
## Primary splits:
## poutcome splits as LLRL, improve=91.42517, (0 missing)
## contact splits as RRL, improve=72.66632, (0 missing)
## month splits as LLRLLLLRLLRR, improve=48.09435, (0 missing)
## housing splits as RL, improve=39.70613, (0 missing)
## pdays < 28.5 to the left, improve=36.09150, (0 missing)
##
## Node number 4: 20854 observations
## predicted class=no expected loss=0.04090342 P(node) =0.7788027
## class counts: 20001 853
## probabilities: 0.959 0.041
##
## Node number 5: 617 observations, complexity param=0.02403846
## predicted class=yes expected loss=0.4246353 P(node) =0.02304216
## class counts: 262 355
## probabilities: 0.425 0.575
## left son=10 (177 obs) right son=11 (440 obs)
## Primary splits:
## duration < 2.5 to the left, improve=44.239210, (0 missing)
## month splits as RRRRLRRRLRRR, improve=13.032810, (0 missing)
## housing splits as RL, improve=10.480080, (0 missing)
## job splits as LLLRRRRRRLR, improve= 5.212144, (0 missing)
## pdays < 199.5 to the right, improve= 5.145055, (0 missing)
## Surrogate splits:
## contact splits as RRL, agree=0.724, adj=0.040, (0 split)
## campaign < 6.5 to the right, agree=0.720, adj=0.023, (0 split)
## pdays < 606 to the right, agree=0.716, adj=0.011, (0 split)
## default splits as RL, agree=0.715, adj=0.006, (0 split)
##
## Node number 6: 5095 observations, complexity param=0.01236264
## predicted class=no expected loss=0.3022571 P(node) =0.1902752
## class counts: 3555 1540
## probabilities: 0.698 0.302
## left son=12 (1500 obs) right son=13 (3595 obs)
## Primary splits:
## contact splits as RRL, improve=54.21724, (0 missing)
## month splits as LLRLLLLRLLRR, improve=36.53465, (0 missing)
## housing splits as RL, improve=28.85552, (0 missing)
## pdays < 371.5 to the left, improve=16.16480, (0 missing)
## job splits as RLLLRRRLRRR, improve=14.23330, (0 missing)
## Surrogate splits:
## month splits as RRRRRRLRLRRR, agree=0.848, adj=0.483, (0 split)
## campaign < 24.5 to the right, agree=0.706, adj=0.002, (0 split)
## balance < -891.5 to the left, agree=0.706, adj=0.001, (0 split)
##
## Node number 7: 211 observations
## predicted class=yes expected loss=0.2227488 P(node) =0.007879897
## class counts: 47 164
## probabilities: 0.223 0.777
##
## Node number 10: 177 observations
## predicted class=no expected loss=0.2768362 P(node) =0.006610151
## class counts: 128 49
## probabilities: 0.723 0.277
##
## Node number 11: 440 observations
## predicted class=yes expected loss=0.3045455 P(node) =0.01643201
## class counts: 134 306
## probabilities: 0.305 0.695
##
## Node number 12: 1500 observations
## predicted class=no expected loss=0.1893333 P(node) =0.05601822
## class counts: 1216 284
## probabilities: 0.811 0.189
##
## Node number 13: 3595 observations, complexity param=0.01236264
## predicted class=no expected loss=0.3493741 P(node) =0.134257
## class counts: 2339 1256
## probabilities: 0.651 0.349
## left son=26 (3347 obs) right son=27 (248 obs)
## Primary splits:
## month splits as LLRLLLRRLLRR, improve=46.61047, (0 missing)
## housing splits as RL, improve=18.72056, (0 missing)
## pdays < 371.5 to the left, improve=12.98456, (0 missing)
## job splits as RLLLRRRRRRR, improve=10.95079, (0 missing)
## balance < 1597.5 to the left, improve=10.86097, (0 missing)
## Surrogate splits:
## day < 1.5 to the right, agree=0.933, adj=0.028, (0 split)
## pdays < 472.5 to the left, agree=0.932, adj=0.008, (0 split)
##
## Node number 26: 3347 observations
## predicted class=no expected loss=0.3274574 P(node) =0.1249953
## class counts: 2251 1096
## probabilities: 0.673 0.327
##
## Node number 27: 248 observations
## predicted class=yes expected loss=0.3548387 P(node) =0.00926168
## class counts: 88 160
## probabilities: 0.355 0.645predicted.model6 <- predict(heart.tree, newdata = test.x, type = "class")
matrix.conf <- table(Class=test.y,Predicted=predicted.model6)
percent.correct <- 100 * sum(diag(matrix.conf)) / sum(matrix.conf)
print(sprintf("L'error de classificació és: %.2f %%",100 - percent.correct))## [1] "L'error de classificació és: 9.15 %"mosaicplot(matrix.conf)
confusionMatrix(predicted.model6,test.y, dnn = c("Prediction"))## Confusion Matrix and Statistics
##
## NA
## Prediction no yes
## no 10153 944
## yes 106 273
##
## Accuracy : 0.9085
## 95% CI : (0.9031, 0.9137)
## No Information Rate : 0.894
## P-Value [Acc > NIR] : 1.29e-07
##
## Kappa : 0.3072
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.9897
## Specificity : 0.2243
## Pos Pred Value : 0.9149
## Neg Pred Value : 0.7203
## Prevalence : 0.8940
## Detection Rate : 0.8847
## Detection Prevalence : 0.9670
## Balanced Accuracy : 0.6070
##
## 'Positive' Class : no
## Obtenim una predicció molt similar, quasi del 0.91 amb un Kappa del 0.3072. La predicció del ‘no’ és molt bona, 0.9149 i em millorat la del yes a un valor de 0.7203
Cart:
fancyRpartPlot(heart.tree, caption = NULL)
Aquest arbre és molt més interpretable que el de l’algoritme C5.0. CART ha realitzat l’eliminació d’atributs i només ens hem quedat amb duration, poutcome, month i contact.
Per finalitzar farem servir un altre mètode per classificar basat en la teoria que vam veure a clustering. El que fem és trobar valors veïns similars a partir de les distàncies dels diferents atributs, és a dir, troba similaritat de posició amb N dimensions. Hem de definir un valor de K. Aquest valor indica el nombre de veïns per a definir un nou grup o en aquest cas predir el resultat.
Preparem les dades:
# Knn
bank.data <- read.csv('bank-clean-discret.csv', sep = ',')
bank.data <- bank.data[,c("poutcome","month","pdays","day","housing","duration","contact","loan","balance","y")]
# separem el 70% de conjunt d'entrenament del 30% per a avaluar
random.data <- read.csv('random-data.csv', sep = ',')
random.data <- unlist(random.data, use.names=FALSE)
train <- bank.data[random.data,]
test <- bank.data[-random.data,]
# en les dades d'entrenament, separem les dades per a generar l'arbre de la dada objectiu (target)
train.x <- train[,1:9]
train.y <- train[,10]
# en les dades de test, separem les dades per a generar l'arbre de la dada objectiu (target)
test.x <- test[,1:9]
test.y <- test[,10]
## creem la funció per normalitzar
ni <-function (x) {(x -min (x)) / (max (x) -min (x))}
# passem a numèric i normalitzem tant test(x,y) com train(x,y)
train.x.norm <- train.x
train.x.norm <- lapply(train.x.norm[1:9],as.integer)
train.x.norm <- as.data.frame(lapply(train.x.norm[1:9],ni))
train.y.norm <- train.y
train.y.norm <- as.numeric(train.y.norm)
train.y.norm[train.y.norm == 1] <- 0
train.y.norm[train.y.norm == 2] <- 1
test.x.norm <- test.x
test.x.norm[1:9] <- lapply(test.x.norm[1:9],as.integer)
test.x.norm[1:9] <- as.data.frame(lapply(test.x.norm[1:9],ni))
test.y.norm <- test.y
test.y.norm <- as.numeric(test.y.norm)
test.y.norm[test.y.norm == 1] <- 0
test.y.norm[test.y.norm == 2] <- 1Calculem un valor optim de k. Veins propers.
i=1 # declaration to initiate for loop
k.optm=1 # declaration to initiate for loop
for (i in 1:75){
predicted.model7 <- knn(train=train.x.norm, test=test.x.norm, cl=train.y.norm, k=i)
k.optm[i] <- 100*sum(predicted.model7 == test.y.norm) / length(predicted.model7)
k=i
}
plot(k.optm, type="b", xlab="K- Value",ylab="Accuracy level") 
Veint la gràfica optem per a fer servir un k entre 20 i 40
predicted.model7 <- knn(train.x.norm,test.x.norm,cl=train.y.norm,k=25)
print(sprintf("La precissió de l'arbre és del: %.2f %%",100*sum(predicted.model7 == test.y.norm) / length(predicted.model7)))## [1] "La precissió de l'arbre és del: 90.59 %"matrix.conf <- table(Class=predicted.model7,Predicted=test.y)
percent.correct <- 100 * sum(diag(matrix.conf)) / sum(matrix.conf)
print(sprintf("L'error de classificació és: %.2f %%",100 - percent.correct))## [1] "L'error de classificació és: 9.41 %"mosaicplot(matrix.conf)
confusionMatrix(predicted.model7,as.factor(test.y.norm), dnn = c("Prediction"))## Confusion Matrix and Statistics
##
## NA
## Prediction 0 1
## 0 10108 929
## 1 151 288
##
## Accuracy : 0.9059
## 95% CI : (0.9004, 0.9112)
## No Information Rate : 0.894
## P-Value [Acc > NIR] : 1.3e-05
##
## Kappa : 0.309
##
## Mcnemar's Test P-Value : < 2e-16
##
## Sensitivity : 0.9853
## Specificity : 0.2366
## Pos Pred Value : 0.9158
## Neg Pred Value : 0.6560
## Prevalence : 0.8940
## Detection Rate : 0.8808
## Detection Prevalence : 0.9617
## Balanced Accuracy : 0.6110
##
## 'Positive' Class : 0
## k-Nearest Neighbors: